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TABLE OF SYMBOLS 


Quantity 
cross sectional area of metal 
cross section area of fluid 
inner radius of pipe wall 
half width of channel (slab geometry) 


b ba" , used in equations (IV-29) 
and (fV23 1} 


ha/k = Fluid Biot number defined 
by equation (V-33) 


¢ha/y |. = Fluid Biot number defined 
by equation (V-62) 


hbjx = pipe wall Biot number 
defined by equation (V-33) 


outer radius of pipe 
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hE} FrCeA, = parameter defined by 
equation (IV-2) 

pet FnCm4m = parameter defined 
by equation (IV-5) 


Fluid specific heat 
pipe wall specific heat 


coefficient defined by equation (IV-57), 
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coefficient defined by equation (IV-67), 
dg=0,, 1 


Units 
ft” 
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ste 


ft 
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ft 
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h(@, n;t) 


L(x) 


K (x) 


heat transfer coefficient at fluid-wall 8 
interface Btu/ “hr F 


heat transfer coefficient defined by 


equation (II-15) Btu/ “pro 
heat transfer coefficient defined by 5 
equation (II-9) Btu/ft*hr°F 
Function defined by equation (IV-45) 

modified Bessel function of first kind 

of order n 

complex integral defined in Appendix A 

equation (5) 

Bessel function of first kind of order n 

pipe wall axial conductivity Btu/ft hr°F 


pipe wall transverse conductivity Btu/ft. hr°F 


pipe wall isotropic conductivity Btu/ft. hr°F 


fluid axial conductivity Btu/ft. hr°F 


fluid transverse conductivity Btu/hr ft OF 


fluid isotropic conductivity Btu /hr OE 


coefficients defined by equation 
(IV-40), i= 0,1,2,3,4 


coefficients defined by equation 
(UN 242) ig 3.4, 5 


modified Bessel function of second 
kind of order n. 


perimeter Et 
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m RX 
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pipe wall temperature 


radially averaged pipe wall 
temperature 


Time 


temperature functions used in 
Section III 


fluid velocity 


temperature function defined by 
equation (III-21) 


temperature function defined by 
equation (III-92) 


axial coordinate 
transverse coordinate 
Bessel function of second kind of 
order n 
GREEK LETTERS 
fluid thermal diffusivity 
pipe wall thermal diffusivity 
tiny (used in Section IV) 
x/V (used in Section IV) 
fluid temperature 
radially averaged fluid temperature 
("%-t)/€ (used in Section IV) 
fluid density 
-t (used in Section IV) 


fluid temperature function at interface 


ft 
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fe oe 
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x - function defined by equation (V-43) 


pipe wall temperature function at 
interface oe 
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Also, additional notations are occasionally introduced for 
temporary use in the analytical developments which follow. In such 
cases, the new notations will be adequately defined in the text. This 
also implies that the same symbol may have different significance in 
different portions of the text. There should be no confusion arising 
from this practice. In particular, different notation is frequently 
introduced for dummy variables of integration. Also, the symbols 

Pw Sn An’ Pn’ etc., used as eigen values, do not appear 
in the foregoing list. 

We also remark that equations are numbered serially in each 
section. In referring to equations, the serial number alone will 
indicate the correspondingly numbered equation in the same section 
as the reference; otherwise the section number will also appear. Thus, 
in Section V, for example, the reference (12) is to equation (V-12) 
and the reference (II-7) is to equation (II-7), that is, to equation 7 of 
Section II. A reference such as IID2 is to the text of subsection IID2 


and not to an equation. 





ABSTRACT 


Transient temperature distribution in axisymmetric systems of 
incompressible fluid flow through semi-infinite, externally insulated 
cylindrical pipes is considered under the assumption of constant uniform 
velocity and constant material properties. The general problem, includ- 
ing all significant heat transfer mechanisms, except radiation, is 
reduced to an integral equation containing only one unknown function, 
namely an interface temperature as a function of time and axial position. 
By neglecting one or more heat transfer mechanisms, several other 
formulations are obtained. These are Reainre so as to provide 
approximate and closed form solutions not previously available. 

Similar formulations and results are obtained for the case of fluid 


between parallel plates. 
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I. INTRODUCTION 


A. PRELIMINARY REMARKS 

The problem of transient temperature fields in systems of fluid 
flow through circular pipes or through channels formed by parallel 
plates is of considerable practical and theoretical importance. On the 
practical side, one may consider, for example, the problem of the 
calculation of thermal stresses in main coolant piping systems of 
nuclear power plants in which abrupt changes in power output may 
cause very rapid variation in the temperature of the coolant fluid 
leaving the reactor. In such cases one may be concerned with the time 
history of the coolant temperature at some critical point, suchasa 
pump or valve, downstream from the reactor. 

From a theoretical point of view, the problems implied by such 
pipe flows are of great interest. They involve two dependent variables, 
the temperatures of the fluid and of the pipe wall, each of which is a 
function of three independent variables: time and two space coordinates. 
These functions must satisfy appropriate partial differential equations 
and external boundary and initial conditions and must also conform to 
each other in an appropriate way at their common interface. Such 
problems are called ''conjugated.'' Few solutions for such problems 
are given inthe literature. Those which were previously known usually 


involve drastic simplifications such as considering only steady state 
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behavior or restricting space variation to one dimension, or considering 
the pipe to have infinite thickness, or neglecting important mechanisms 
of heat transfer, or a combination of these. 

Accordingly it is important to consider such problems inas gen- 
eralized a formas possible both for the purpose of obtaining compact 
alternate formulations and also for the purpose of obtaining closed 
form and approximate solutions to problems in which the degree of 
simplification is less than that involved in previously known solutions. 

It is believed that the work reported herein is successful in both 


regards. 


B. STATEMENT OF THE GENERAL PROBLEM 

It is appropriate at this point to make an explicit statement of the 
general problem which will be considered herein. We are concerned 
with the flow of an incompressible fluid axially down the length ofa 
semi-infinite cylindrical pipe whose cross section is a circular annulus 
of constant finite wall thickness. It is presumed that axisymmetrical 
conditions prevail so that the space variables are radius r and axial 
position x from the accessible end of the pipe. The outside of the 
pipe is presumed to be perfectly insulated so that there is no transfer 
of heat through this surface. (This condition is approximately satisfied 
in many cases of technical importance. ) 


It is assumed that the physical properties of both the fluid and the 


pipe wall are constant -- that is, they do not vary with respect to time, 
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temperature, or position. We presume that the thermal conductivity 
of the materials involved is orthotropic; that is, we provide for 
different values of conductivity in the radial and axial directions. This 
permits later specialization in different applications, either to the case 
of isotropic conductivity or to the cases where one or more conducting 
mechanisms are suppressed. Subject to this introduction of orthotropic 
conductivity, it is assumed that Fourier's law of heat conduction holds 
in each medium. In most of the work done herein, we provide for an 
interface surface heat transfer coefficient to account for surface 
resistance or other effects. This coefficient may be taken as infinite 
if there is no such resistance. However, inall cases it is assumed to 
be a constant which does not depend on time, position, or temperature. 

We have assumed the nature of the fluid flow, namely, we have 
assumed that the velocity of each element of fluid is the same, a con- 
stant velocity V in the axial direction x. This is frequently called 
slug flow. It is recognized that this is a particular case which is at 
variance with the flow conditions which would actually obtain under 
any circumstances. However, it is not far removed from the flow 
conditions which would prevail with fully turbulent flow and it does 
provide a necessary simplification which has permitted us to obtain 
useful results. 

We have confined attention to one case of temperature variation 
which can be described as follows. It is presumed that the fluid has 


been entering the pipe at zero (or reference) temperature for a 
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solution for the case of cylindrical geometry. In fact, in some cases, 
when certain simplifying assumptions and procedures are introduced, 


the two cases may become one and the same. 


Geo URVEY OF LITERATURE 

The study of heat transfer problems in systems of fluid flow 
through pipes dates back to 1883 when one of the most celebrated 
problems of heat transfer was solved by L. Graetz [11].* In Graetz's 
problem the velocity field is taken to be known, and heat transfer 
takes place under steady state conditions with pipe wall temperature 
distribution specified a priori. In recent years other problems of 
fluid flow have been treated by many authors. These studies include 
such effects as heat transfer by axial conduction and mechanical heat 
generation due to friction, i.e., viscous dissipation. There is an 
abundance of literature in this area of heat transfer engineering. 
However, these problems are concerned solely with behavior of the 
fluid; the action of the containment is considered only in the specifica- 
tion of boundary conditions for the fluid. 

However, the transient case of temperature field in systems of 
flow through pipes or ducts when the wall thermal state is not specified 
a priori, and it is required that its effect on the temperature field in 


the fluid be accounted for, clearly fallsin the class of conjugated 


*Numbers in square brackets refer to the Bibliography, page 154 
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problems of heat transfer. These conjugated problems are more 
involved inasmuch as the mutual thermal effects of the wall and the 
fluid must be accounted for. Accordingly their treatment presents 
serious mathematical difficulties. Even when the velocity field is 
assumed to be known, the thermal problems are still of great difficulty, 
mainly because of the coupling that is introduced by the condition at 

the common boundary of the fluid and the wall. 

At present there are numerous publications dealing with conjugated 
problems of heat transfer. Inall cases, however, they seem to be 
particularized and as yet no general method has been obtained for 
treatment of these problems. 

An early treatment of conjugate problems of heat transfer was 
Pen by T. L. Perelman [25] who considered steady state heat transfer 
to a fluid flowing around a body containing internal heat sources. He 
also considered the asymptotic solution to the type of integral equation 
that occurs in the analysis of such conjugated problems: namely, a 
Fredholm equation of | second kind with singular kernel and infinite 
range. He considered slip flow around a body and also considered 
the laminar boundary layer flow over a thin plate. Perelman neglected 
the axial heat conduction in the fluid. 

N. A. Avdonin [4] considered transient behavior in a similar case 
but with no heat sources present. He also neglected the axial conduc- 
tion in the fluid. He solved the problem by using Green's functions 


and Laplace transforms. He treated the conjugated problem by 
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reducing it to a Volterra integral equation which he solved using 
Laplace transforms. 

D. E. James and N. G. William [17] considered the steady state 
problems of Poiseuille -~- Couette flow between parallel planes. In 
this work the boundary planes were taken to be of finite length in the 
direction of flow and of finite thickness perpendicular to the ite 
of flow. These authors assumed an arbitrary temperature distribution 
at the fluid-wall interface and solved the energy equations separately 
for each material for a constant interface temperature distribution. 
Then they used Duhamel's superposition integral to obtain the general 
solutions interms of the arbitrary interface temperature distribution. 
Next applying the continuity condition for the heat flux they obtained an 
integro-differential equation for the interface temperature distribution. 
They solved this by the method of collocation. 

Recently A. V. Luikov, V. A. Aleksashenko, and A. A. 
Aleksashenko [22], considered conjugated steady state heat transfer 
in circular and plane tubes for the case of developed Poiseuille velocity 
distribution allowing for mechanical energy dissipation. They state 
that general methods have been developed for both steady and unsteady 
heat transfer, but they give no example or discussion for other than 
steady state cases. They essentially reduce the problem to a singular 
integral equation for the unknown temperature at the fluid-wall 


interface. This is a very interesting paper. 
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E. M. Sparrow and F. N. DeFarias [33] investigated the laminar 
flow through parallel plate channels under the assumption of sinu- 
soidally varying fluid inlet temperature. This work is in the spirit of 
the new results to be reported herein inasmuch as a solution for 
sinusoidally varying end temperature can be related to that for step 
change in end temperature and conversely. Accordingly it is important 
to delineate clearly what simplifying assumptions they employed. They 
neglected axial conduction in both fluid and pipe and they assume pipe 
temperature varies with axial position but not with transverse position. 
They considered constant fluid velocity as we do herein. It happens 
that these simplifying assumptions do not coincide with any set 
treated in this thesis. 

V. S. Arpaci [3] considered a transient problem of the type in 
which we are interested. Assuming as we do, the case of slug flow, 
he uses the term "rod'' to describe what we call the fluid. He shifts 
viewpoint in an obvious and nonessential way by considering the pipe 
to be moving at constant velocity with its contents at zero velocity. 
Axial conductivity is neglected in both substances. Radial conductivity 
femeccounted for in the fluid. Radial conductivity of the pipe wall 
is taken as infinite. The transient effect is due to sudden and uniform 
heat generation in the contents. Like ours, his pipe has perfect 
external insulation. He used the method of finite Hankel transforms 


to reduce the problem to the solution of an integro-differential equation 
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which was solved by Laplace transforms. For large values of time an 
asymptotic behavior of the temperature function was found, again using 
Laplace transforms. 

Other significant contributions to the study of conjugated heat 
transfer problems have been made by M. D. Kelleher and K. T. Yang 
[19], R. Siegel [34], and R. Siegel and J. M. Savino [35]. 

F. E. Moreland [23] employed numerical inversion of Laplace 
transformation to deal with heat exchanger problems which are similar 
to our problems, except for suppression Of all heat diffusion inthe fluid. 
D. D. Hiep [13] extended this work significantly by including the effect 
of axial heat diffusion. (His treatment is for porous media, but our 
formulation could be similarly cenemmaeed by use of some fictitious 
material and geometrical properties.) His study contains a large 
bibliography. The problem considered in Section VI hereof is a special 
case of that considered by Hiep. We are able to obtain analytical 
results whereas Hiep's analysis is purely numerical. 

Two earlier solutions which provide significant points of departure 
should also be mentioned, even though idealizing assumptions made 
the analyses very simple. One of these is by W. Munk [24] who 
assumes identical constant temperature in both pipe and contents at 
any cross section. He neglects all diffusion of heat and arrives ata 
solution which repeats the time-temperature history at the entrance 
point, but with a propagation velocity less than the velocity of the 


fluid. The second, to which many investigators have made contributions, 
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is discussed in section 15.3 of [7] where many references are cited. 
Arpaci [2], page 353 et seq., gives a very compact treatment of this 
problem in which there is no heat diffusion in either the axial or radial 


direction. 


D. SUMMARY AND OUTLINE OF SUCCEEDING SECTIONS 

The following telegraphically brief statements are intended to 
provide a summary of the content of the succeeding sections. 

Section II - Derives basic equations which will be employed in 
subsequent sections. 

Section III ~- Considers radial and axial heat diffusion in each 
material. Arrives at integral equation (III-70) for slab geometry and 
integral equation (III-181) for cylindrical geometry. No solution of 
these equations is attempted. 

Section IV - Radially averaged temperatures are considered in 
each material. An ordinary perturbation solution is constructed 
which is valid "behind the (temperature) front.'' A singular perturba- 
tion solution is developed which gives variation ahead of the front. A 
scheme for matching these at the front is given. 

Section V - By averaging the fluid temperature in the transverse 
direction and neglecting axial conduction in the pipe wall,closed form 
solutions are obtained in the case of slab geometry: equations (V-39, 
40). The case of cylindrical geometry turns out to differ only as 


regards the definition of a certain auxiliary function so that the 


solution has the same expression as for the case of slab geometry. 
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Section VI - By taking radial averages for both fluid and solid 
temperatures and by neglecting axial conduction in the solid, a 
simplified version of the preceding problem is obtained. (This may 
also be obtained by neglecting axial conduction in the solid in the case 
given in section IV). A closed form solution is given: equation (VI-21) 
emacs ( VI-23). 

Section VII - Here the problem is greatly simplified by assuming 
(in effect) that there is no radial variation in temperature. Axial 
conduction in both materials is considered, however; (i.e., we take 
average fluid temperature equal to average solid temperature at any 
section and employ the equations of section IV). The resulting 
solution, (VII-10) is not new but is included here for completeness and 
comparison purposes. 

Section VIII - This section includes discussion and recommendation 


for further study. 
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Il. FUNDAMENTAL EQUATIONS 


A. GENERAL REMARKS 

All the analysis presented in this thesis relates to a definite 
physical problem which will be described shortly. The real physical 
problem is so complicated as to preclude the possibility of dealing 
with it without making certain idealizations. In what follows immediately, 
the problem will be described in mathematical terms and the idealiza- 
tions which are common to all the subsequent analysis herein will be 
delineated. The simplified problem which results is still of great 
complexity. It will be further discussed in section III. No solution of 
this problem will be presented but it will be shown that the problem 
can be cast ina more compact manner which gives some promise of 
being amenable to analysis. 

In Dies to obtain additional results, it has been found necessary 
to make additional idealizations or simplifications. These lead to 
several simplified problems which are specifically delineated and 
dealt with in sections IV through VII. For all of this work, however, 
the derivation given immediately below and the equations which result 
therefrom are fundamental. 

1. Physical Model 

The basic physical problem of interest is that of heat transfer 


taking place between a conduit and the fluid flowing therein. The former 


Be 





is perfectly insulated on its exterior surface. It may be either a 
circular pipe with uniform wall thickness or a duct composed of two 
parallel plates. (In what follows these two cases will be referred to 
as "cylindrical geometry" ee "slab geometry, '' respectively. It is 
assumed that the fluid entering the conduit has had constant temperature 
for sufficient time that the conduit and its fluid contents have reached 
this constant temperature throughout. Then the temperature of the 
fluid entering the conduit is suddenly changed by a finite amount and 
is thereafter maintained at the new temperature. (The end surfaces 
of the conduit also experience the same step change in temperature. ) 
It is required to predict the unsteady temperature fields in both the 
fluid and the conduit for all locations at any subsequent time. 

Figure 1 illustrates the two geometries which are considered. 
That for the cylindrical geometry is self-explanatory, except to remark 
that there is no variation whatsoever with respect to angular coordinate. 
The case of slab geometry involves identical parallel plates which are 
semi-infinite in extent in the x direction, which is the direction of 
fluid flow,and are infinite in the z direction. It is assumed that there 
is no variation whatsoever with respect to the coordinate z of any 
quantity. In making a mathematical formulation of the problem, it 
will be convenient to take advantage of the symmetry with respect to 
the plane y = 0; for this reason the auxiliary coordinates Vy and Me 


are used as indicated. 
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[eux = OS 


a - cylindrical geometry 


70 X=o0 


b - slab geometry 





Figure |; Geometry and Dimensions of Conduits 





2. Assumptions 

a. The fluid is considered to be incompressible. 

b. The fluid is assumed to be flowing with constant velocity V 
(slug flow). 

c. Material properties are taken to be constant and independ- 
ent of time, position, and temperature. 

d. It is assumed that there may be an additional thermal 
resistance (possibly due toa very thin layer of foreign material at the 
interface) the heat transfer coefficient of which is a constant. 

e. Viscous dissipation is assumed negligible. 

f. Heat transfer by radiation is considered to be negligible. 

g. Fourier's law of heat conduction is taken to be valid 
throughout. 

h. The diffusivity appearing in Fourier's law for fluid is 
assumed to include not only the molecular diffusivity but also an eddy 
diffusivity which is effective. The combined diffusivity is assumed 
to satisfy assumption c. 

i. There is variation only with respect to time, the axial 
dimension, anda radial dimension, Figure la,or a transverse 


direction, Figure lb. 
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B. DIFFERENTIAL EQUATION 
Pee Ovlandrical Geometry. 

There is no difficulty in establishing, by customary methods 
of derivation, the following partial differential equations: (1) which 
describes the fluid and (2) which describes the solid pipe wall. We 
have anticipated simplifications which will be made in later sections 
hereof by considering artificial media in which thermal conductivities 
in the axial and transverse directions may differ. These conductivities 
are distinguished by a second subscript, a to denote axial or t to denote 
transverse. If we are considering isotropic materials (as we doin 
section III hereof, for example) we take the axial conductivity to be 
equal to the transverse conductivity. However, if in order to make 
mathematical headway ina certain analysis we neglect axial conductivity 
in the solid (say), we can obtain the corresponding mathematical des- 
cription simply by setting K = Q. 

ma 
There should be no difficulty in recognizing the significance 


of the various terms which appear in equations (1) and (2) below. 


08 ee 
kpo 38 - eM ax * Kees (3 +)" FS 3 ia 
OT 1 2 (.2T OT 
K Stk, to (rd)- 5) AC er, (4) 
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One may observe that (2) can be obtained from (1) by taking 
V =0, replacing 98 by T, and replacing subscript f by subscript m. 
2. Slab Geometry 
The partial differential equations governing the temperature 
distribution in the flowing fluid and in the solid duct wall are 
Ky 22 + Ky, Dye Fe vee a 2 i 
X - ¢ ot 


for the fluid, and 


ma 3, OT 4k mt a 
2 


for the pipe wall. Observe the use of subscripts 1 and 2 for the 


coordinate y, consistent with figure lb. 


C. INITIAL AND BOUNDARY CONDITIONS 
1. Cylindrical Geometry 
aa fOr ek haitc 
(1) Initial Condition, t = 0. The initial temperature of 
the moving fluid is considered to be uniform, and without any loss of 


generality, it may be considered as zero. 


LIN 


8@=0 0£x < © ; O Qa . t=a0 


WA 


ae 


ITN 


(2) Boundary Condition at x =0. At this boundary the 


entering fluid temperature is given a sudden step change to a higher 
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value thereafter. Without loss of generality, the amount of increase 


may be taken to be unity. 
Q=1 RO 5) 10'S GO eee 


(3) Boundary Conditions as x—»@. As x approaches 


infinity the fluid temperature decreases to its initial value. 


0 —>o X—»0o, orga , t20 
(usually, this condition is employed in the form: @ remains bounded 
as X—>oo.) 
(4) Boundary Condition at r=Q. It is required thatas 
r— >» 0 the fluid temperature remain finite. This condition may also 
be stated by requiring that the temperature distribution in the fluid 


to be symmetric about the center line. 


Fither 
ee @(v,.x,t) < 00 o<x < 00 5) t20 
r—-o = 
or 
08 = r=0 » ° “<x ¢< @ 2 20 
or 3 


b. For the Pipe Wall 
(1) Initial Condition, t = 0. The initial temperature of 
the pipe wall is considered to be uniform and equal to the fluid tem- 


perature at the initial time. 


yT=o ofxdm , asr gb 5» t=O 
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(2) Boundary Condition at x = 0. The pipe wall is 


assumed to experience the same unit step change of temperature as 


that of the fluid. 


a 
T 
i) 
ad 
Tf] 
@) 
i) 
tA 
KH 
IN 
oOo 
cr 
ny 
@) 


(3) Boundary Condition as x ~»oo. Similarly it is 


required that the pipe wall temperature approach its initial temperature 


as x approaches infinity. 
T—o x-—*o , a<7r Qo Sts (8 


(4) Boundary Condition at r = b. Since it is assumed 
that the pipe wall is perfectly insulated at the outside surface, then 


the Fourier heat conduction law requires that 


Ol =o 0o£x <0 ; Sol |e it 60 
or 


c. For the Interface Between Fluid and Wall (r = a) 


Continuity of heat flux is required at the interface between 


fluid and pipe wall. 


Contemplating a constant heat transfer coefficient h at 


this boundary, we may thus write 


“Ky: 08 =h(e-T) 

or Oeste GG ft - a to 
ON era g = 1) 
mt ar 
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It is well to note at this point that if the heat transfer coefficient h 


approaches infinity, then the above conditions reduce to equality of 


temperature and heat flux, in the following manner. 


0 =F 
0<¢X <m , =a 
K 09 . OT 
ft or mt or 


2. Slab Geometry 


The details and physical significance of the choice of initial 


and boundary conditions for this case are similar to those discussed 


for the case of cylindrical geometry. We will therefore state them 


mathematically as follows. 


a. For the Fluid 


(1) 


(2) 


(3) 


(4) 


imptial Condition £ = 0, 


O= 0 OLX Ko, ogy sa 
Boundary Condition at x = 0. 
6=1 X=0 geo SUG 


Boundary Condition as x—»oo. 


Q=0 


x—? 2 SPE ss 


Boundary Condition at yy, =4- 


9 O$x ¢o > y= a 


ral 
c©o iD 
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ier. figure lb; the symmetry abourthe face y, =2 provides this 
condition. ) 
b. For the Wall 
(1) Initial Condition, t = 0. 


T =0 o<x <0 » OS Y> 


T 
oF 
cr 
| 
O 


(2) Boundary Condition at x = 0. 
hs I x = 0 ; 04 y2<b ; t>o 
(3) Boundary Condition as x — oo. 


T =0 x —» 00 ; 0<¢y,Eb , t>o 


(4) Boundary Condition at y> = b. 


we 
O42 


ec. For the Interface Between Fluid and Wall Vq = Yous 


o<x <o ; Yo =b , tz2o 


t 
oO 


eee 


mt 
84. 


h (6 - T) 


OSx<0O | ¥) -¥2-0.1=9 
eOee 1) 


i] 


Observe that the algebraic signs given here are a con- 
sequence of the fact that V1 and Y> are opposite in direction; Cf. 
Figure lb. Again it is easy to observe that as h approaches infinity 


the interface condition above may be written as 


oO 
WA 


x €O 5, Y,=¥n =O, tZzo 


31 





: 
- 
= 





D. SPECIALIZATION OF ENERGY EQUATIONS 

In what follows particularizations of equations (1) and (2) for the 
cylindrical geometry and equations (3) and (4) for the slab geometry 
are introduced, which involve various techniques for simplifying the 
energy equations. These methods involve radially averaging the fluid 
and the conduit wall temperatures by integrating with respect to radius 
(for cylindrical geometry) and with respect to the transverse co- 
ordinates (for slab geometry). Thus appropriate definitions of average 
temperatures and surface film heat transfer coefficient are introduced. 
Further simplifying assumptions are made, suchas neglecting the 
axial conduction in the conduit wall; these simplifications give rise to 
interesting problems. Also since the difference between cylindrical 
and slab geometry disappears when both the fluid and the wall tem- 
peratures are averaged, only the cylindrical geometry is considered 
in this case. 

1. Cylindrical Geometry 

a. Radial Averaging of Temperatures 


We consider isotropic conductivity in both media: 


Kea Sa Ke (>) 
Kaan nt = Ke (6) 


We also define the average temperatures 
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7,8 


. [rede a 
gn- 2. 2 feeds C7.) 


a 2 
a 
[ede . 


0 


and 
b 


[eter b 
=, eaaeioueas d d 
T ie fet 5 (8) 
ror a 
a 


al 





and we introduce an "equivalent" surface film coefficient 


x 20 ol 
p = - Ke (SF Jose = ai Gaae 
aaa t aan se 


Note that h* relates heat flux to difference between average 
fluid temperature and average wall temperature. 

Then, by integrating equations (1) and (2) with respect to 
radius, between appropriate limits, and introducing the definitions 


above, we find 


2 x 
6 vm lt 20" _ (e.T*) (10) 
ox? a ox me ot AK, 


for the fluid mean temperature, and 


vr" aT 


Oxcueawerot 





A’P (o% T") (11) 


m Km 


po 
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for the pipe wall mean temperature. Equations (10) and (11) are 
essentially those considered by Hiep [13]; although his equations 
appear to provide for more generality (in order to apply to porous 
media), a change of variable shows that they ese mathematically the 
same. The introduction of h” may legitimately be objected to on the 
grounds that physically it would be quite difficult to evaluate. However, 
its use provides a degree of generalization compared to the physically a 
reasonable simplifying assumption that radial variation of temperature 
in the fluid is so small that e* = 6 (r=a)= T(r=a). Moreover, 
even in cases where h-—~ oo, the use of h” <0 permits taking into 
account some resistance to radial heat flow in the fluid. 

The initial and boundary conditions for (10) and (11) follow 
naturally from those of (1) and (2). 


(1) Initial Condition t = 0. 


Bo 5 PP sg 0 <x <oo : t=o0 
(2) Boundary Condition at x = 0. 

6 aT a1 x =0 ; t>o 
(3) Boundary Condition as x — oo. 

ee 510 : T—»+o x—00 : t>o 


b. Case where h”™ Approaches Infinity 
An interesting special case may be obtained if it is assumed 
that radial diffusion is so great that the temperature is essentially 
constant throughout any particular cross section. This is accomplished 


by letting h* approach infinity. Then from (11) we obtain 
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Keer take oT oe 
(e*-7") = > (oe Sued (12) 


Je 


Substituting for (@ - T*) from (12) above in (10), and noting that as h”™ 


approaches infinity, then @ = Ie, yields 


x (1+ Kae ) So" y BE a (peta ok (13) 
t K, Ag / ox? J. A, S, Ot 

for the common fluid and conduit wall temperature distribution. The 
initial and boundary conditions for (13) follow naturally from those of 
(10) and (11). 


(1) Initial Condition, t = 0. 


i 
>) 


aie = 0 o<x < 00 ’ t 
(2) Boundary Condition at x = 0. 


T=) x=0 , ie 


ny 
© 


(3) Boundary Condition as xX 49. 


T*—»r 0 x—>oo , t 


ny 
O 


The solution of this system is treated in section VII. 
c. Case where K, = 0 
By neglecting axial conduction in the solid pipe wall, a 


tractable system is obtained. We rewrite equation (10) with no change 


2 Ps x 
ve _Y 26 ~ 1 oe eloag 10 





Multiplying equation (11) by Kym and allowing Ke to approach zero, 


wae id S 





aan’ _ 


The initial and boundary conditions are 
(1) Initial Condition at t = 0. 
acon ee 0fXKm 4, t=0 


(2) Boundary Condition at x = 0. 


al 
“< 


oe | X¥=0 2 t>a 
(3) Boundary Condition as x —»go. 

—_ O XK —P 00 oy t>o 

The solution of this system is dealt with in section VI. 
d. Another Special Case 
If radial averaging is done only for the fluid and variation 

with respect to radius is retained in the solid pipe wall an interesting 
problem results. However, this situation leads to an integro- 
differential equation the importance of which does not seem to be 
commensurate with its difficulty. Accordingly an additional simplifica - 
tion is made: axial conduction is neglected in the solid pipe wall. To 
accomplish this we set K,., = 0 in equation (2) and we also write 
Kmt = Ky. For the fluid we write Kr, = Kg = Ke. Note that the 
radial effect is absorbed in the radial averaging for the fluid. It is 
necessary to define a aa equivalent heat transfer coefficient pes 


follows. 
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-k, 27 . -K, - £(0" T) at rea (15) 


or 


Then equations (1) and (2) with 9” still defined by (7) reduce to 


% 
ee OPAL } a 
ie 
x “f Xx X “ Ap K. 


and 


af us (17) 
fi an 
(17) is exactly the form taken by (2) with Kya = 0. (16) is the same 
as 10) except for the new definition of surface film coefficient. The 
equivalent heat transfer coefficient hi , defined by (15) relates inter- 
face heat flux to the difference between conduit surface temperature 
and fluid mean temperature; compare with h” defined by (9). Clearly 
oo ~~ 
h <h . The reasons for introducing h are covered in the remarks 
following (11). 

The appropriate boundary and initial conditions are: 

(1) Initial Condition. The initial conditions for equations 


(16) and (17) are exactly the same as those for equations (1) and (2) 


o* = 0 O¢K (00 » t=0 
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(2) Boundary Condition at x = 0. The same conditions as 
stated for equations (1) still apply to equation (16): 
psu X=O 4» to 
We observe that since equation (17) is explicitly 
andependent of x, it is therefore not required to specify any condition 
at x = 0 for equation (17). 
(3) Boundary Condition as x-»g@. It is again required 
that @* remain finite; viz. 
O- = 0 KX —» 00 ; t >o 


(4) Boundary Condition at r=b. Assuming perfect 


insulation at the outer surface of the cylinder, then 


or 
Tyan yob 3 t 


Nv 
° 


(5) Boundary Condition at r=a. At the interface of 

fluid and pipe wall solid the continuity of heat flux requires 

~K,, 28 wh (ot T) r=4 9 tZe 

m = 
or 

We note here that if h approaches infinity, then 
the condition of continuity of heat flux above reduces to equality of 
temperature at the fluid and wall interface. 

This system is analyzed in VB1l. See discussion at 


end of IID2 for explanation of why the system of IID2 is treated first. 
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2. Slab Geometry 
For all the simplified cases in Section IID1 in which average 
temperatures in both fluid and pipe wall are considered, the geometry 


of the situation is fully described by the parameters A AL and P. 


f? 
Thus, specializing the case of slab geometry in the same ways leads, 
as one may easily verify, to exactly the same equations. However, 
the case treated in JIDld still retains radial variation in the pipe wall. 
Accordingly, we expect that the slab geometry counterpart will be 
described by slightly different equations, which we now proceed to 
obtain. 

a. Differential Equation 


Treating equations (3) and (4) for the case of flat wall 


duct similarly we define 
a 


a 
rh thos (18) 
J d 


Then considering isotropic conductivity in the fluid, 


integration of (3) yields 


z ¥ x £ 
SY Mot S +(e 7) (19) 
Ox 4 * Ox K. Oot aK. 


for the fluid. Similarly with the assumption of zero axial conductivity 


in the wall, we obtain 


a9 





2 
Coast | (20) 
doy Z <m ot 


rp 


for the solid wall. Note the introduction of h in (19); see remarks 
following (17). 
b. Fluid; Initial and Boundary Conditions 
(1) Initial Condition t = 0. 
9* = QO O¢x <0 ; t = 0 
(2) Boundary Conditions at x = 0. 


ae 1 x = 0 : t 


hy 
S 


(3) Boundary Conditions as x—>»oo. 


sie 


6 —> 0 x —» OO E t > 0 
c. Solid Wall; Initial and Boundary Conditions. 
(1) Initial Condition, t = 0. 
T =0 Os$ Y2 <b ; t = 0 
(2) Boundary Condition at Y2 =b. 
OT =" y, =b ; t >0 
O92, : : 
(3) Boundary Condition at ¥2 = 0. 
-K_0T = % (6 - T) v= 0 9 tz 0 


Oy> 
This system which is similar to that obtained in subsection IIDIld is 


somewhat easier to deal with. Accordingly, it is treated first, in VA 


while that for the system IIDId is dealt with in VB. 
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Ill, INTEGRAL EQUATION BORN Eater 


PaeelNT RODUCTION AND PLAN OF ANALYSIS 


Equations (II-3) and (II-4), with Ke, = Ky = Kp and Kya = K =K 


mt m 


f 


so as to provide isotropic thermal conductivity, involve two dependent 
variables, 6(x,y,t) and T(x,y,t), each involving three independent 
variables, along with a formidable number of boundary and initial 
conditions. Similar remarks hold for the case of cylindrical geometry. 
These are complicated systems to deal with and we have achieved 
solutions which are useful for calculating purposes only by making 
various simplifying assumptions. These solutions are ee in other 
sections of this thesis. Inthis section we will represent the system 
of equations (II-3) and (II-4) in another form which involves only one 
unknown function of two independent variables. No additional conditions 
are required for the complete specification of this unknown function. 
The form of this equation, (70), is that of an integral equation of what 
is knownas "first type'' (that is, the unknown function does not appear 
outside an integral sign) with complicated kernels and forcing functions. 
Although the method of procedure is essentially the same for the 
case of cylindrical geometry as for the case of slab geometry, the 
details are sufficiently different to justify treating the two cases 
separately. We treat the case of slab geometry first and in subsection 
IITD we return to the case of cylindrical geometry. In each case the 


method of procedure is essentially as described in the next paragraph. 


4] 





We define an unknown function } (x,t) which represents the 
(unknown) temperature of the fluid at the interface and a similar 
function y (x,t) which represents the unknown temperature of the 
solid at the interface. For the fluid we solve the problem of finding 
Bite interior temperature 0(x,y,t) when the surface temperature is 
given by P(x, t), subject to other appropriate boundary and initial 
conditions. It is convenient to consider this solution as the sum x 
two solutions, one involving a homogeneous boundary condition at the 
interface and accommodating the boundary condition at the end, x = 0, 
and a second solution involving a homogeneous boundary condition 
fone — 0 and the unknown function P(x, t) at the interface. The 
solution to the latter problem is expressed in terms of an integral 
involving the function @(x, t); it is obtained by use of a procedure 
which is due to A. N. Lowan [20,21]. We proceed in the same way 
to express the solution for the solid temperature T(x, y,t) in terms 
of the unknown surface temperature y (x, t). The interface heat 
balance condition provides two relations, one of which may be used 
to express the function W (x, t) in terms of the function ¢ (x,t). 

The other relation then becomes the integral equation for the unknown 


function p (oot is 


B. AUXILIARY PROBLEMS (SLAB GEOMETRY) 
We will deal with the problem of the temperature distribution 
in the fluid first. Then, by taking the flow velocity as zero and by 


using subscript m (for the solid) rather than f (for the fluid), and 
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by changing the designation of the surface temperature distribution 
from ¢ (x,t) to y (x,t), we can very quickly obtain the corresponding 
solution for the temperature distribution in the solid. We immediately 
write 
O(x,y,t) = uj(x,y,t) + uo(x, y, t) (1) 
where us(x,y, t) are the palerions to the two auxiliary problems 
considered in the nae of this subsection. 
1. Fluid; Homogeneous Boundary Condition at Interface 
For peonincity in this subsection we will omit subscripts, 

and will replace them when we combine solutions: & will denote 

Ke, u will denote uj, y will denote y, (Cf. Figure lb). Thus the 


problem of interest is 


z 2 
o4 ,ou _VY uw lt du (2) 
Oy xe *K dx *K at 
The boundary and initial conditions are 

u=1 x =0 O<y4a t> 0 (3) 
u-~0 x -» OO O<y<a t2 0 (4) 
u = 0 0<x < y= 0 t>0 (5) 
OuU_ oO 04x <0 yrza t>0 (6) 

oy = ~ 
eae 0<x <0 O<ys a t= 40 (7) 
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we obtain 


Using an overbar to denote Laplace transform with respect to time 


Ba 8a Vv wt 
oy? Ox? 4 


Oo 
pad 


Ss CL 
or 
This incorporates the initial condition 
which becomes 


All the boundary conditions 
will be the same in Laplace transform form except that at x = 0 
U= 


Thus assuming a solution of the form 
U=zAe ~ 


sin( y) 


7) 
we observe that all conditions are satisfied except that for x = 0 
and the condition of zero slopeat y=a 
the characteristic values 


Imposing the latter yields 
fn = (n+1/2) TT /a 


Because of linearity we can write 


(10) 


00 
aa A,X ; 
u(x, 4 5} = » A, e 


Ain ( p. 4) 
n2=0 


Substitution of (11) into (8) yields 


An = 





+ V(¥) re 


xX ——» 00 


(12) 
The negative root is discarded in order to assure boundedness as 


of trignometric functions yields 


The condition at x = 0 and use of the orthogonality property 
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A, = “S é 
nt (13) 


Substituting from (13) and (12) into (11) we obtain 


= — ¥2 
U(x54,5) \ Ba in (Hy) <i E -x \IP t Zz (14) 





where 
z 2 Zz 
V 
= “+ aa 
Mt oe) (15) 
The inversion of (14) may be obtained using standard tables 
[29]. Thus 
0 
U(x, y, t) = > Yj (y) X (x, t) (16) 
where eee 


Y, (4) = — (17) 


1 ele XP pit 
aie ) = tele. erfe (SE ie ee serge (Stet i“ (18) 


O 2yxt 
In these expressions, although there is explicit dependence upon the 
integer n, this is not reflected in the notation onthe left. This is 
done merely for simplicity of notation and will also be done with 
similar definitions to follow. In every case the symbols are used 


under a sign of summation with respect to n. 
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For future purposes we evaluate the partial derivative of 


(16) with respect to y and evaluate at y = 0. 


| le: YXoue) (19) 


nN=0 


where 
/ 
Y (O)=1 (20) 
} 
As the reader may come to appreciate, when the many preliminary 
results are assembled into the final equation, a compact notation is 
convenient. Functions Y, and xy and constants Ha are the first of 
several such notations which will be introduced. 
2. Fluid; Non-homogeneous Boundary Condition at Interface 
In what follows Ustands for Uj (x,y; t) and it is governed by 
the same equation (2). However, the following changes are made in 
the boundary conditions (3) and (5) 
= 0 (rather than u = 1) x = 0 O<yga t>0 
U= P (x, t) (rather than u = Q) 0<x<00 y=0 t20 
where | > (x, t) | £M<o@ forallxandt. We further make the 


following substitution in order to drop out the convection term. 


2 
u(x, y,t) = W(x, y,t) exp [Ys is ve | (21) 


Then equation (2) and the appropriate boundary and initial conditions 


reduce to the following equation and conditions for W(x,y,t). 
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Ox? OY? ~ «& Ot (22) 
w=0 X=0 4 ofysa, t>o (23) 
W < 920 X—rag » OLY CA , C20 (24) 
W = 9 (x,t) o<x<o , 162 2 t>o (25) 
OW _ o<xKo0 > Yaa y tro (26) 
OY 
W= 0 afX€co » O84 ca , t=0 (27) 

where 
vx , vt 
= -~VAa VE 28 
(x,t) = (x,t) exp pw + a (28) 


Again taking Laplace transformations with respect to t and accounting 


for the initial condition we obtain 


SN, Raw . sie 
ox? oy? x en 


All the boundary conditions will be exactly the same as above except 


(25), which becomes W = g (x,s). 
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Following the method of solution used by Lowan [20,21] we 


consider the function 
oo 0 
W (x,y, 5) = # J [Foosa.s) nent hin ue dude | (30) 


where 


F (o,y;)) = cosk [\(a-y)] (31) 


Cos (Aa) 


and 
2 2 
= — 6/2 
(oo (32) 
Note that ha and iE are dummy variables of integration. Direct 
substitution shows that (29) is satisfied. It is easy to verify that all 
the boundary and initial conditions are satisfied. In particular we 


note that 
@ 00 
W (x, o,s)= 2 ff (5,5) Mm ps Pon bX dpa = 4 (x,S) 
970 
by the Fourier identity. As regards boundedness 
|W(x2998)] $ 2 fin px Aim pe dp dt = a < 00 


This completes the proof that (30) is the desired solution of (29) 


which satisfies all the auxiliary conditions. 
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Continuing to follow the procedure described by Lowan [20, 21], 


we define 


ferEcasy.t) de = Flosy; ») = FO.y sf S ) (33) 


re... FE, (a,y;t) is the inverse of F (a,y;A). The subscript P 
indicates the dependence upon the parameter BP; the dependence of F- 


upon s is made clear by expressing A in terms of s, as above. The 


convolution theorem then gives 


00 00 1 
W (x54, t)= Ef [rin px rinpy dp. de fate ysVglsstendan (34) 


Further defining 


it 
G.(a,4,5¢) = fe (2545 t)dt (35) 
then 
OG (a,y,t) - 
eee Seg 7S ; 
as be (% 3 t) (36) 


Thus substituting from (36) in (34) yields 
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@0 00 Gs 
Wxge) = Ff fain pr acnpe dudg (2 (82) 9¢¢ t-n)dy (37 
4 J en 


If we now take the Laplace transform of (35), we obtain 





£ | % (4. 456)) = Ld [fi (ay;t) = = F(4,4;a) (38) 
= i Yt) 
S Z(s) 


where 


Yis)= cosh [a(a-y)] 5 Bis) = cosR (MA) (39) 


The singularities of (38) are simple poles corresponding to 
the zeros of cosh( i a); of course, these are functions of the variable 
_ ye 2 4 4 : 
p . Note that cosh (Aa) = 1 + a/2!+ Na*/4! +... which converges 
everywhere in the finite plane. Thus, although ) itself involves a 
fa | 
Square root , Yd does not, and there is no difficulty to be en- 
countered dealing with branches. Thus, the residue theorem affords 


the desired inverse of (3 8) 


00 5.6 
Gu (a,y;t)-1O, mice (40) 


Z (0) | 4 
n=O 
0s 5= Sa 
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The summation in (40) extends over the roots of 
cosh (Aa) = 0 (41) 


namely the numbers 


»n - pis 8 (n+ yz) SH oo py, eA Ippeee. OL) 


(Do not confuse the constants fn of (10) with the continuous variable 
P of (30)). Thus 
z 2 . 
ae =~ [(m- 1/2) tr + P") 
ae 
(43) 
2 2 
zs OM (Ps + ) 
where S, are the poles of (38) in excess of the simple pole at S = 0. 


Note that they are functions of the variable — Therefore the formal 


inverse of (38) is 


-« (pap )t 
Gul y3t) = costlv(a-a)]_ 2 nyt "(n+ M4) cos [oa(o-y)]} € = (44) 
CosR (Ma) a! ee eee 


The following derivative is needed in (37) 


(psp yt 
0 Gp(a.y;¢) = - ye -1). Pn Cos [Hx (2-3) ] € (m j (2) 
ot N=0 


pul 





Hence, substituting from (45) in (37) we obtain 
aye p esfrncoalff ey 9) 
(46) 


W (x.y, t) = 


jer" Dun po % vo we dp [andy 


O 
In order to complete the analysis, substitute (28) in (46) and then 


Thus we obtain 


substitute the result into (21) 
(ve 


eae at) oS se) fn cos | | Pa (a aif [8O° 2 y) Anu 
(47) 


CO 2 

-P 
ane a ie a 
so that we 


We are able to carry out the integration with respect to a 


finally obtain 


S.2) oo 
U(x,4,t) = » yc ff PCE, t-1) Kilex, tn jende (48) 
N=0 fd 


where 
(49) 


Vlad = Z (-1)" cos [Pa (o-y)] 
. (V4 op?) 





TE Pe Be (X- A 
K, (xe stale ie ot & 8 y 
(50) 
, i fs | Se 
4a 4a 


ft 


Se 





For later purposes we need | ot which is evaluated below. 
gz0 


es 
an fas. t-r) Ki (xu ¢e,t,n) dude 2) 


In order to assure that there is no misinterpretation, we 
remark that the integration with respect to € is from zero to infinity 
and that with respect to n extends from zero tot. 

3. Combination of Solutions for the Fluid 


Simple addition of (16) and (48) yields 


oo oo t 
O(x, yt) ‘) (ene 5 yotan f 4G) K (studs (52) 


Likewise by superimposing (19) and (51) we obtain, 


00 


5 |, LI botd ef [#052 Kyler tonanas | (53) 


where, having restored the appropriate subscripts, 


Weg) = =" (54) 
ae 4 “iy per 

~ | *“F X-2X%-H, rec ( Xt2X% pst as 

A (xt) = Le Be “erge( PRES) 4 ¢ -erfe( Xt20 a (55) 

Y, 43.) —s = [Pata (56) 


216. 





Also, of course, (10) and (15) define My and ro 





(57) 





It is well to pause at this point to review the physical nature 


of the problem for which (52) provides the solution. 


below. The differential equation is 
26 v2 se _ lr 
ax% ae ‘Ox dy oT t 
The auxiliary conditions are 
6 = 1 X= 0 5) 0 £1 {a 
6 < a0 X —> 00 3 Des es 
@= $(x,t) ofX <a , y= 0 
08 | é 7 
ai Ose SECTS) g,* 
O=0 O0¢ X <0 J o$g,¢a 


4. Solution for the Solid Wall 


This is given 


(23) 


t20 

t>o 
t2o 
tro 
t=o 


The solution for the similar problem for the solid wall is 


easily obtained from the preceding solution for the fluid. 


All that is 


necessary to dois to write T in place of 6, Y in place of ¢ , replace 


subscript f by subscript 


m, replace a by b, take V = 0, and 


introduce a new symbol A. (defined below) in place of #,,. 
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Oo 0 a) 
T (x1g,t)= » 0 Xz (x,t) + Yow frye, t-1)K2(x.8,t7)2148 | (58) 
n=O es 


where 
Y;, (42) = Gamcrh i : (59) 
~dnX* AnX 
Py a X-2% Xt i x+2% Ant 
Xz (x,t) " e erge( er jre cnet) (60) 


vey os [An(4 -y) 
Bi, = (1 cos Lan ( jai 


2 
z “XA [ 2 2 
An Se | =(t- -(€+ 
Kote VF ee re vp tel) on Ed (62) 


An = (nei) 2 (63) 


(Do not confuse d\n above with X , used temporarily in (42)) 


a0 wot 
gis *) [x 1Osot) af fv gt-n Rute tin) dacs | (64) 
H=0 : 
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Cc. DERIVATION OF INTEGRAL EQUATION (SLAB GEOMETRY) 
Equations (52) and (58) thereby give the fluid and pipe wall 
temperature distributions respectively. However these solutions are 
in terms of the arbitrary interface temperature functions P (x, t) 
and Y (x,t). We now return to the actual boundary conditions that 


prevail at the interface, viz. 


Ke h($¢-¥) Yio 9 ox Seo, tBe (65) 
and 
Knog BOER) diz J ox So > t>o (66) 
2 


(Note: The signs on the left hand side of (65) and (66) are consistent 
with the choice of coordinates used in Figure lb). Solving for yw from 


(65) in terms of P , we obtain 


v= ¢- Se (38) (67) 


as 


Substitution from (53) in (67) yields : 


(exe) ot 
VY (xt)= Pit) - a Y [x09 +f P(E ,t-) (68) 


n=0 


Ky 00 §>t, 1) 4ads) 


56 








Furthermore, substitution of (68) in (66) yields 


00 ot 

re 

«(2 =k, [X(t + (Et) K(x, sendy) Te 
oe OY. q=0 if j I] ‘S i 


nMz0 


Substitution of P (x,t) from (68) into (64) and the result into (69) 


yields 
00 ot = 
Kk) 2 (x5 +) of [Ka lasgstsn)- $(¥,t-7) - ay (5, t-n} 
N=o oa ete 
a) t-y 
+f foe, t-n-f ) K(E,«, t-1 2 4p de dade + 
oO oe 
e0 ot 
kK) XxX, (x,t) +f fect K (xe. tr) dy4e = 0 
“or a) 


It should be clear that we could, with equal facility, express the 
problem in terms of My (x, t) inate of ? (x,t) as above. The 
distinction is trivial. 

Equation (70) is the desired integral equation, the only unknown 
in which is the function } (x,t) which represents the surface tem- 
perature of the fluid. 

Equation (70) is considerably simplified if there is no surface 
resistance. Letting h, the heat transfer coefficient, approach 


infinity, reduces (70) to 
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Op co t 
Gi= \ et) a ra X,(*>t) + | faeeeu- 


Neo 


Ki(og, ton) Me Kaley, tn) ]dady (71) 
f 


fomen = Oo . 


D. AUXILIARY PROBLEMS (CYLINDRICAL GEOMETRY) 

For the case of cylindrical geometry, the treatment is quite 
similar to that which has preceded. However, it is not possible to 
derive the equations for the pipe wall by simple manipulation of the 
equation for the fluid as was done above. Thus it is appropriate to 
consider four auxiliary problems similar to those dealt with before, 
two for the fluid, arriving at equation (165), and two for the solid, 
arriving at equation (170). These equations will then finally be 
combined so as to arrive at the desired integral equation. By suitably 
redefining the several functions which are involved, it will be possible 
to express this equation, (181) exactly in the form of (70). 

1. Fluid; Homogeneous Boundary Condition at Interface 


Consider the problem mathematically stated as 


Bu_vUw, 1 qu, 4 _ 1 w 72 
Ox? &K x YY @ drt X Ot site 
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with the boundary and initial conditions: 


r= 1 X=0 3 o <r ga 3 tro (73) 
U=0 X—> 0 » °§ rS$q 5 tro (74) 
nyee od re) a 

= ‘ , s  tze (75) 
Us o OSK <0 9 r=H g t>o (76) 
U= o 0K <coo lg OK SA, ge (TT) 


Assuming a solution of the form 


U = AJ.(é@r) X,(« >t) (78) 


and substituting (78) in (72) yields 


2 2 
Oe NN. Oxi = | OX (79 
ox Xx dx PX aoe , 


Equation (78) already satisfies condition (75). Condition (76) yields 
the characteristic equation 

J,( 6%) = 0 (80) 
the roots of which are denoted by Bo n=1,2,... . Because of the 


linearity we can write 
Oo 
us \ And, (far) Xi (x2) (81) 
n=} 


The solution of (79) is obtained by taking Laplace transform with 
respect to time, which yields a second order homogeneous ordinary 


differential equation the solution of which is 


a9 





X (xt) = = B oe [Ye a Gar. ae] (82) 


In (82) the exponential with positive exponent has been discarded so 
that condition (74) may be satisfied. Substitution of (82) in (81) and 


absorbing B into Ay , we obtain 


Us ) Aad Oan) <x E ~x\(%) +R) + 3 (83) 
n=4 


Thus at this point we observe that the inverse of (83) satisfies (72), 
(74), (75), (76) and (77). Condition (73) provides a means of evaluating 


Ay. Taking x = 0 (73) gives 


Co 
= E ) An Je (Por) (84) 


n=4 


multiplying by rJ, ( f.Y )dr and using the orthogonality relation of 


Bessel functions; 
ao 
JI( Pat )To(Par)rdr = 2 == oT ( (An) (85) 


if n=m and zero if n#m, one easily finds 


l 


eee 
J,(B.%) 





ae (86) 
p, 


Thus we now have 


doy Ti (Br ee 5 or |e VT Pa =| = 


) 2% 
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The inversion of (87) is tabulated [29]; viz. 


oo 


MGs.) = > Rr) Xi (xt) (88) 


na 


where 


Racer) = 1 Fe (Bar) 


(89) 
Pn J, (6%) 
men. iF [o" “t : 
Xy = = K+ 204,% 
) Pe erfe ( H2iee ra eae verge (22 ) 90) 


For future purposes we need the partial derivative of (88) 


with respect to radius evaluated at r=a. 





Ze 


Fluid; Non-homogeneous B.C. at Interface 

The problem considered here is the same as that in the 
preceding subsection except that the following changes are made to 
the boundary conditions (73) and (76). 


Ul 


i 


0 (rather than u = 1) X=0 


uU= 9 (x,t) (rather than u = 0) ape A r=& » tro 


where \¢ Bee t) | <M < 00 for allx,t. A substitution of the form 


2 
Caemxat je avy (it x,t) exp (Ye - am 


i” Ge (92) 
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reduces (72) to 


2 
OW , 1 OW OW - | ow (93) 
and all the boundary and initial conditions will remain the same 
except that at r =a we have 
z 
EE aS 
= —— + = 
W(a,x,t) = P (x, t) exp ( 5K GK g(x, t) (2) 
Taking Laplace transform of (93), we obtain 
= as ee 
Oi |\WwW,wWwW@iosy 
eee a OS 


(95) 
Again all the boundary conditions will be the same except (94) which 


in Laplace transforms will become 


pe es (96) 


Following the method of solution used by Lowan we consider 
the function 


where 


W(r,x,5) = 2/ [Fay 9 (E15) Aim MX Men HE dpdt (97) 





F A, ane = I, (Ar) 
( ) (aa) (98) 
and 
ve p+ >. (99) 
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Note that # and € are dummy variables of integration. Direct 
substitution shows that (95) is satisfied as well as all the boundary 


conditions. In particular once again we note that 


Wlanslee ff 5(e,5) Aiea Ux Rim WE Ap dE 


(100) 
= 9 (x5) 
by the Fourier identity. As regards boundedness, 
oo 00 
NV < Zt c ° _ 2M 
\W(r,%,5)| ¢ 2H I F acm pr Mon WE dpe = 22 ¢ a (101) 
00 


This completes the proof that (97) is a solution of (95) which satisfies 
all auxiliary conditions. Let us now define 


Qo 


Est ia Z - 
€ bi (Qr;t t= F(a, 5d) = F (a4; {ps ) (102) 


oO 


i.e. , ae (a, r;t) is the inverse of F (a,r; ). Note again that the 
subscript ie indicates the dependence upon the parameter Ne = ee 
dependence of F upon s is clearly indicated by expressing X in 


terms of s, as shown in (102). The convolution theorem then gives 


0 0 3 
Wirsxyt)= 2 f frinpr fin PE op de ic on)g(e.t-n)ér (103) 
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Further, defining 


i 
Gp lasrit) = fFp(asrst)dt poe 
O 
gives 
0Gp (ar; t) = Fe (a,r3t) (205) 


ot 


Substituting from (105) in (103) yields 
t 
W (onthe af fain dMnne dpde Sela) alg -tndan (106) 
oO o 


If we now take the Laplace transform of (104), we obtain 


a Gp Carst)]- Enc e;t) dt 


(e) 


5 ot [hE (ars6)] (107) 


5 
ie = F(a, er; d) 
ae: | ace 
5 2(s) 
where 
Y(s) = I, (0d4) d Z(s)= I, (ra) (108) 


64 





Using the (temporary) substitution 4 = Lp , (98) may be written 


= J(2r) 
J, (8) 


(109) 


However Jo(3 ) is an entire function having only real zeros. Thus 
the only singularities of F are poles at the real roots, fin of 
Jo( Ba) = 0 (110) 


The corresponding values of s are: 


ous _ (p+ p: ) (111) 


: 


Gp asrst) 


Following Lowan, we write the desired inverse as 
|, Sn 


cE 
3 aS oy 
a (H+ Bn) 


I.(pe) Bade Irn) e 
ar (pa) tau (Cree?) FT (Pra) J ( Px) 


Thus taking the derivative of (112) and - into (106) yields. 


Wltsxst) = 4.1 4, ZolBar) ott a (F,t-2)- 
ka 7 J (Ana) J 


= (113) 
: er ae dnde 


e) 


(112) 


It is possible to carry out the integration with respect to p 


in (113) and obtain 
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a o t 


Uirwtd= \Rucedff ey .tn) Boopt nares cs 


n={ 0° 
where 
J, (Ant) 
Rip Ge) (115) 
Ba, (Poe) 
x po (x-¥) a(%e«th) 
Ky (ny tin)= (=o ean . e 
1 CEEG) 


? Zi 
Exp Susy) -x) - oP ~(g +x) 
4xN) ax vy 
Note that Kx, Boe ee Inge ae ele B Oy ecl it Bewve mietyorya Ky(x,¥ PE) 1 hee ae 
by (50) only in the eigenvalues; otherwise they are identical. 
For later purposes we need the partial derivative of (114) 


with respect to r, evaluated at r = a, 
OU. aaa 
Se] =) ffeG-tuKiey tase cin 
eee 


3. Solid Pipe; Homogeneous Boundary Condition at Interface 


Consider the problem stated mathematically as 


e 


(118) 


(Pe) = 27 


Y 
ct 


ox? Y or 


with the boundary and initial conditions: 
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U= 1 X=o , 464 (6b > = ae (119) 
“=o X—> co 3 ALrSb > beNe (120) 
ou eo O<KX Loos r=} » tro (EZ) 
te= 5 0< X<00 6 f= a » tro (122) 
Wo oSxcoo » alr fb , tee (123) 


We first define the following notational device for a linear combination 
of Bessel functions of the first and second kind of an integer order. 


Let 


cn (rb: ¥) = J.( se) ¥(¥b)-Y, (vr) J, (8b) (124) 
Nowassume a solution of the form 


Y=) Al 6, (163 5) ea x, t ) (125) 


Substitution of (125) in (118) yields 


Dy Z OX 
— eae _t om (126) 


The function (125) already satisfies (121), i.e., 
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Pete ke) = ,6,(b,b;%)} =0 Gia 
or r=b 


Condition (122) yields the characteristic equation 
06, (2%, 65%) =0 (128) 


the roots of which are denoted by YY, , n=1,2, ... . Because of 


linearity we can therefore write 
| Oo 
uo = ) Aa 6 (12 63 b) X2(x,t) (129) 
n=] 


To solve (126) we take Laplace transform with respect to 
time. This yields a second order ordinary homogeneous differential 


equation, the solution of which is 


Nol S) = 6B <&xp [=x ' wae $1] (130) 


In (130) the exponential with positive exponent has been discarded so 
that condition (120) may be satisfied. Substitution of (130) in (129) 


and absorbing B into AY we obtain 
00 
“(r,x,s)= ) An 6,(+,b;¥) = KY Yq + 5/K 131 
ay ) n es ag am exp Unt ( ) 
n=| 


Thus at this point we observe that the inverse of (131) satisfies (118), 
(120), (121), (122) and (123). Condition (119) provides a means of 


evaluating A. Taking x = 0 (119) gives 


68 





y ie S (r, be Mer) (132) 


It is not too difficult to show that (124) has the usual orthogonality 
properties of Bessel functions with the weighting function r. Thus 


we merely present the result as 


° ae & (a, ey) 
J rBlorivs:6 (eta) © Be EO yyy 
2 a’ 6 (a A, b; Sm) - L Erbe 1) 


ifn = mand zeroifn#m. Use of this orthogonality property yields 


af Gi (a,b; xn) 


a 
> | a E* (a,b: %n)- BC *(b.b:%s) 


(134) 


Substitution of “4 in (131) yields 


© 16,(a,b5 5.) & (Gbs ¥) 
_— ne OA BSB) Er AG D5 BH) | ly? $s (135) 
4 -) a ts (0,6; 8) - b* E> (bb; x) me rex b+ = | 


Note that the parameter s appears only inthe last term. The inverse 


of (135) is tabulated. Hence, 


oo 





Me secs. t) =) R3( X, (x, t) (136) 
where rat 
Ry(r) = & 61 (4b be) 66 we, (137) 
" HE @bibn)- BE (bb; 5) 
- nx A Sn X xt 
Xi(x.t) =e rfc (Apes +2 erge (a (138) 
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The partial derivative of (136) with respect to r and evaluated 


eter = ais 


= / 
Ea = \ 83 (a). Xz (xt) (139) 
or 


r=Q2 a 


where 


L/a 
E Elbe) 
¢ 06 (a,b; ¥n) 


/ 
R3 (a) = (140) 
This result will also be used later. 
4. Solid Wall; Non-homogeneous Boundary Condition at Interface 
We consider the immediately preceding problem with the 
following changes in boundary conditions 


Uu=0 (rather than u 


U 
| ane 
ra 
Tl 
oO 


agrsb tro 


2 


i Y (x,t) (rather than u 


I 
oO 
rat 

tI 


where 





Y (xX, | €M <oo forall x andt. After Laplace transforma - 


tion we have 


ya an 
ou / 2 OY yess oS ae 
Se ee — (pe ] = —-U 141) 
ayer ar | i = ) 
with the conditions 
“=o X =0 acr db (142) 
U< co X —» 00 acrcéb (143) 
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| 


GY 
<7 
QO 


O¢K Ca De ee (144) 


oe 
% 


W(x,5) o<x < @ 5) c= 6 (145) 


The solution of (141) may be constructed following the procedure used 


by Lowan. Thus, consider the function 


00 00 
U = 2 [LF (a,b. 5r) H(G,5) Apr An py dp dy (146) 
00 


where 
(a,b, wjr)= Ted Kb) + Ko) TO6) aay 
I, (Aa) ky(Ab) + Kola) T,(Ab) 
and 
ire pot S/o (148) 


(IT and K, are modified Bessel functions of first and second kind of 

order n). Observe that p and ie are dummy variables of integration. 
Direct substitution shows that (141) is satisfied. It is also 

easy to verify that all the boundary conditions are satisfied. In 


particular we note that 


i (0.2.5) = 2/19 (5,5) Rinp hin rg dp dy aay 


= V(x, s) 


ga! 





by the Fourier identity. As regards boundedness, 


@ oo 
lu Chesik 2 [rin px Ain py dude = an < 00 (150) 
00 


This completes the proof that (146) is the desired solution of (141) and 
satisfies all the auxiliary conditions. Again we follow the procedure 
indicated by Lowan, and the sequence of analysis is precisely that 
employed in subsection IIIB2 for the case of slab geometry; the present 
easceis daiferent only in details. Accordingly, we abbreviate the 


treatment and write immediately, 


se t 
u (r,x, t) Dad fonpx sinppdude | Mier y ¥ (¥en | Geer) 
Oo 


re) 


and 
é ee" | “fa ys) 
ap (a, bor; t) 7 E F(a, betas) |= ae E is | (152) 
where 
Y(s) =] (ar) Ky (ab) + K Orr) J, (44) (153) 
Z(5) = J, (da) K(Ab) + Ko (da) J, (Ab) (154) 
Using the temporary substitution =i¥ , (147) may be written 


E(a,,0;¥) = J, (80) (8b) - Ji(Sb) Yo(Ee) _ ob (tb3x)_ 


(155) 
J, (8a) Y, (¥b) -J, (3b) Y, (8a) 6, ( &, (@.b33) 


(2 





where we have used the definition (124). It is reasonable to assume 
that the function , (a,b; ¥) considered as a function of ¥ has 
only real zeros. Thus the only singularities of F are poles at the 


real roots, Y,, of 


56; (a,b on ro (156) 
The corresponding values of s are. 
2 2 
Sn a-« (¥, +p) (157) 


Following the procedure used by Lowan, for the inversion of (152), 


we obtain 


OO Se 
Cn Siok See Y(sn) € 


— 2 oe (158) 
Z(o) s| 22 | 
ie SS Bee 
Gula, b r;t) = I,(#r) Ki(Hb) + Kol Br) Th (Pb) 
tra) K,(pb) + Ke (pa) I, (pa) ae 


(Xn ee)t 
2) 260(¥,b; 4) € Fee 
En tee [b iS (4,b;%,) + 0 ‘S (a, b5¥) 


Thus taking the derivative of (159) and waite into (151) yields. 


u(nx,t) = ) Rat fre. t-n)Ka(xseston) Inde = (160) 


where 


i6 








Rau(v) = _ Vn& (vr, b; a) 


A NEN Sci 2920 a (161) 
bE, (2b; ¥n)+ aE (a,b; a) 


2 
-AXTn 2 


z 
K GE tn)= t= Exp P| - exp| ae] (162) 





For future use we require the following 


co t 


0 / 
ou] Z 
ae) Ralf fv 1) Kilx.g, to rerde (163) 
where 
‘ x" 
; : 96,(4, b; ¥n) ~_ 
5 
6, (a, b; Sn) 


5. Combination of Solutions (Cylindrical Geometry) 
aeeeetior the Fluid 
Simple superposition of the solutions (88) and (114) yields 


00 ot 
BG xt)= | RUe) Xi(x.4) + RCO [ete euKies, oars (165) 


n=} 


where having restored the appropriate subscripts 
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iii) =Ri(r) = a Jo ( Bx) 
Pn J, (par) 





(166) 


yx Rex 
beis.t)= |. et | -¢ X= Dn Ket e . er Xt DX te 
Se Se did eax: jee”. i retrare rx 67) 


\ 


Ms (x-¥) = (%g Ba + V /d4o%) 
f 
2 - € 


K,(x.4,t.1) = \%e 


(168) 


rite 


2 2 
= -xX £ 
exP (§-x) — <KP (§+x) 
40% ¢ 4x6" 
Similarly superposing (91) and (117) we obtain 


= x, (Ce) ar Peete Kk, (x,€5 ton)éeds| ite) 
| oD ff 


Gee bor the Solid Pipe 


Again simple superposition of (136) and (160) yields 


O90 ot 
T(e,x,t) =: ) Rse)Xte t) + Role) {Fy Kilgore (170) 
N=; o 0 
where again restoring the appropriate subscripts 


Rae) = /6,(4b3¥n),6, (¥,b; ¥n) est 
Bn 0" 6°@,b3%)- B6,"(b,b 5¥u) 


Ra(r) =-=& (ry bs hm) (172) 
yi be, (a,b; ¥n)+@ © (A,b;%n) 
2 
2 -% 0" 
- [Xm _om € 173 
K(x g.t4) = T b © re (17 ) 


oat) —_ 
(n-ne 


Me 





similarly superposing (139) and (163) yields 


00 aot 
if 
| - | = Y [Rao X (xt) + Rie Fycpten)kibonta es (174) 


= 


where 


/ { 
3, Cease (175) 
7 ( 2 06, (b, b; Sn) 4 
0. 5 (a, b; Sn) 


Z 
¥n (i776) 
b 6, (a, b; Sn) 
6, (a,b; Sn) 


\ 


f 
Ry (2) 


+ 8 





BE. DERIVATION OF INTEGRAL EQUATION (CYLINDRICAL GEOMETRY) 
We now proceed as we did in section IIIC for the case of slab 

geometry, Equations (165) and (170) give the fluid and solid temperature 

distribution each in terms of the arbitrarily preassigned function 


(x,t) and W(x,t). If we now return to the actual boundary condition, 


namely that of continuity of heat flux at the interface, we have 


m0 P. - | 
“KL | &(?-¥) (177) 
Yaa 
and 
ol) = 
BE] fh (p-y) (178) 
reQ 
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and solving (177) for y in terms of + , we obtain. 


ve at 22 | + > (179} 
£ OF es 


Substitution of (169) into (179) yields. 


= > - “| X,(xot) + if P(Est-n) K(x, Suis (180) 


nal 


Further substitution fory from (180) above into (174) and the result 


into (178) yields 


K ce Ry (a) Xp (xt) + (0) ffatoss >”) 
aetn)- £5 fe pen effete): 


m= 


co a 
“Ky(€.6, ty, x) oe pe KY Deloot) +f [05 .t-a)K lug, ts1) reso 


n={ 


(181) 


The same remarks as those following equation (70) for the integral 
equation of slab geometry may be made here too. Similarly therefore 


if we let h approach infinity, we obtain 


20 oo t 
yk "= Ry(a) Xp (ust) + Xy( (ot) ffecee) 


n= (182) 


[KiQog, bute R, (Canes ern tn) isp ne 
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This is equivalent to (71). Note that (182) can be written more 


compactly as 


> [x “othe [foreaK' eS otrrlendg |= (183) 


N=( 


where 


X (ut)= ue R, CON Ce ea) (184) 
£ 


¥ / 
K (o8st.7)= Kj (x's t2) 4 C Ry (2) K,005,61) (185) 


A similar abbreviation may be made for (71). 

The reader should note that (70) and (181) are of exactly the 
same form; their only difference is in the particular form assumed by 
the kernel and driving function and in the fact that the counterparts of 


R3(a) and Ry (a) appearing here are simply unity in the case of (70). 


F. REMARKS ON SOLVING THE INTEGRAL EQUATION 

It is obvious that the integral equation (70) or (71) for the case of 
slab geometry, or (181) or (182) for the case of cylindrical geometry 
will not be solved easily. It involves multiple integrals, one with 
variable upper limit and one with fixed infinite upper limit. The kernel 


and driving functions, although well behaved, are quite complicated. 
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One of the ways of dealing with such a complicated problem is to do 
what has been done inthis thesis, namely, make additional assumptions 
so as to simplify the problem to the point that it becomes analytically 
tractable. Our attention has been exclusively directed in this way and 
we have made no serious attempt to deal with these fundamental integral 
equations. However, the following remarks and conjectures may be of 
interest. 

All the analysis herein shows that while ¢$ (x,t) is a smooth 
function, almost all its variation takes place in the immediate 
neighborhood of the wave front x = Vt. It is conjectured that a contour 
map of > (x,t) would be similar to that shown in Figure 2. In this 


figure the contours are shown as straight lines merely because the 


x=Vt 
a 





Figure 2. Conjectured Contours of © (x,t). 


fie, 





nature and degree of curvature are unknown. It is tempting to consider 
Peplacing the boundary and initial conditions $= 1 along x = 0 and 

? = 0 along t = 0, by the conditions $ = 1 along OA and 6 = 0 along 
OB so as to provide a field in which the "action" of (x,t) is more 
evenly distributed. This might perhaps be useful in an attempt at 
solution via Gallerkin's (or a similar) method in which a solution is 
assumed which involves unknown scalar constants and these are 
determined so as to minimize an error functional. In sucha procedure, 
it is evident that simplifications would result from assuming the con- 


tour lines to be straight. 
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IV. METHOD OF PERTURBATION 


In this section we employ the method of perturbation to treat the 
problem described by equations (II-10) and (II-11) and the associated 
boundary and initial conditions. The reader will recall that these 
equations resulted from having taken radial (or transverse) averages 
of temperatures in the fluid and in the pipe wall. Thus radial (or 
transverse) conductivity is, in effect, absorbed into the surface heat 
transfer coefficient h* defined in equation (II-9). The influence of 
axial conductivity appears explicitly in the following analysis. 

An outline of what will be accomplished in this section may be 
helpful to the reader. First, an ordinary perturbation expansion is 
made interms of the parameter € = 1/V. It turns out that only even 
powers of € are involved. The zeroth order solutions are identical 
to well known solutions described in section 15-3 of [7]. The effect 
of longitudinal conductivity does not appear in this zero order solution 
but it does appear in subsequent solutions. It is clear from the sequence 
of higher order perturbations that the solutions fail at and ahead of 
the wave front, x =Vt. Severe singularities are encountered at the 
wave front and the solutions are identically zero ahead of the wave 
Front. 

Thus a singular perturbation development is introduced to explain 


the behavior ahead of the front. The lowest order solution for this 


Sl 





singular perturbation expansion is similar to the known solution for a 
moving fluid, with diffusion and with an isothermal boundary condition. 
Terms in this singular perturbation development involve unknown 
functions of time which are determined by a process of matching with 
the solutions valid behind the front. 

It should be remarked that the system of equations dealt with in 
this section is obtainable from either slab geometry or cylindrical 
geometry. The process of radial or transverse averaging of tem- 


peratures erases the distinction between the two geometries. 


A. ORDINARY PERTURBATION ANALYSIS 

Equation (II-10) and (II-11) together with associated conditions 
are repeated here for reference purposes. The notation has been 
slightly changed; we omit the asterisks which indicated radial averages 
and we introduce the convenient constants bb and b, defined below. 


In the fluid, we have 


z 
26 106 _ 00 
. & b, (0-7) (1) 


where 


ao 


b= (2) 
€ A, 


+1 0 


subject to the conditions 
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6=0 o< X< 00 2 t=o 


For the solid we have 


Xx oT = ol . ‘le at} 
ot 


where 


subject to the conditions. 


RecA: xX =0 > t20 = 
To X —» 00 » t>o 
T=0 0 <X S00 9 t=0 


It is convenient to normalize with the following change of 


variables 
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(3) 


(5) 


(6) 





See 


N=X/y =€x ; e=l/y 


é 


so that equation (1) and conditions (3) become 


222 28 _ 208 


€ dye an me tb (8-T) 
O=1 (=o  » tz 
O —» 0 \—r0 » tro 
@=0 Of < t=0 


and similarly equation (4) and conditions (6) become 


2 Onibe Ol lt 
KO. 6 2) = SS a 
ie Cis (6 T) 


T=H1 = 0 tz 0 
Wessis "| —» 00 T>o 
T=0 05" <2 t=o 


LOE, 


Equations (8) and (10) with conditions (9) and (11) can be approached by 


the use of Laplace Transformations as was done by Hiep [13], but 


there seems to be no way of solving the transformed equations or 


performing the inversion other than by strictly numerical procedures. 
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This system of equations may be considered as depending on the 
parameter ae which appears in (8) and (10). For large values of V, 
this parameter is small so that it is reasonable to consider a per- 
turbation expansion in € 

(Note that € itself does not appear in (8) through (11), only Pe 
thus the expansion we develop is in powers of a that is, only even 
powers of € are involved. There are other ways of showing that odd 
powers should not appear in this perturbation expansion, but the one 


given here is sufficient.) Thus, 
2 
O(n. t;é)=O (nat) € O(n st)+--- (FZ) 


Ws @)) 2 si je) ee Gist soe (13) 


2 
Substitution of (12) and (13) in (8) and (10) and equating powers of € , 


the following sequence of equations are obtained 





06, . 00 Ve " 

an et * ont 

DO6c . 20 20 

Oe Ora | One i) =a lla 4 16 
ay Ot (eee) = 9 ra —— 


for the fluid and 
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le _ 4, (0,-%) <0 (17) 





ot 
2 
a Sb, (ei-m) = *, Se (18) 
ut 
ah Patt 
Olea lorena |) rc i> el 
= be ( i) Ss One 2 2 (19) 


for the solid. The boundary and initial conditions are 


O(n .t)= 4 ae 2 t>o 
(20) 
8,(n,t) =0 of <coo, to 
. (ns t) = so 2 #20) | 
see (el 
62 (nt) = 0 0S | <o tzo 
and 
To (4 >t) ={ (= 0 J t20 
(22) 
ibs) Se oS<o , t =o 
c(4,t) =o (=o 4 t>o 
t>d (23) 
Tela ne)) ee oSn <oo , t=o 


Since the order of these equations is less than that of the original 
equations, only one boundary condition can be imposed. It is not 
expected that the solutions will be uniformly valid. Starting with (14) 
and (17), we take the Laplace transformation with respect tot. The 


transformed equations are 
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JES Opel, (3,216) ec (24) 
3 


(s+ bz) T = b2 @, (25) 


Expressing T,in (25) interms of 6, , substituting in (24), and solving 


the first order homogeneous differential equation yields 








= b bibz 
5) S(S+ bz) Sed | 1 | “ay 
6, (128) = = =u [on (5+ Pls a ) es 


The procedure for inversion of (26) and (27) is as indicated in 


Example 7-5 of [2]. The results are 


Omit) =o for t<4 (28) 


@, (n,t) = nant belo "ET (ag (ay )ay for t>n (29) 


and 


T, (yt) = 0 for t<n (30) 


Use t)= b. e ye I, (z¥ay) dy for t>Y (Bla) 
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mnere a = by bon and T=t-r. Note the use ofa inthis sense; it 
should not be confused with a linear dimension. It may easily be 
shown that 6, (X,t) and To(%,t) both tend to unity as t-woao. 

It is clear that the two solutions obtained for Tp(\,t) and 6,(4,t) 
have a wave like character with no disturbance for %)t. It is well 
to pause a moment to refer back to the physics of the problem. For 
fixed n (i. e., fixed value of x) no disturbance takes place until t Dae 
or for fixed t, there is no disturbance for the region y» >t. In what 
follows the condition T = t ~-n = 0 is referred to as the ''wave front. '' 
The pipe wall temperature T,(n,t) is continuous at the wave front 
with discontinuous slope, while the fluid temperature 9,(N,t) is 
itself discontinuous at the ae front. 

fie tirst corrections 6,(",t) and T, (y,t) are governed by 
equations similar to those for 6, ( n,t) and T,(n,t) except that these 
new equations are nonhomogeneous with the second derivative of 

G,( yH,t) and Tg(n ,t) as forcing functions. These non-homogeneous 
terms vanish ahead of the wave front, and thus no correction ahead of 
the wave front is expected in the sequence of solutions 05, To, 6;, 
Ti, --- . Inthe next subsection we will demonstrate this by actually 


constructing the solutions 0, and T, and will see that the present 


1 


perturbation expansion fails at the front. 
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B. HIGHER ORDER SOLUTIONS 
The corrections 6,(n,t) and ly (n ,t) are also obtained by 


using Laplace transformation. From (18) 


= 5 = 
(s+ br) Th = bz A + x 9 Te (32) 
i 


and from (26) we have 
2 


OTe _ be b, b2 _ _ by bz 
“ee rere ial 2 exp| n (s+ bi bibs. (33) 








Likewise from (15) 


— Li 


8 + (seb) 8 = bit + 4 





(34) 


and from (27) 





¢ b 
+. (S+ egies . xp |-n (S+bi- he ) (35) 

oye S S+b2 S$+b2 

Substituting (33) in (32) and solving for T in terms of 8, and then 

substituting the result in (34) results ina first order non-homogeneous 


differential equation for 0, (‘> s). 





25 Ova oe lomlali Bn ea hia! 
or tir oe) 4 eel +5) | n(s+by Sry Ce) 


The solution of (36) is easy to obtain. 
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rn y nb bz b w b, bz 
4 ee | ( Sr id Ee u(s+ "1 S+b2 ) a 


It is now merely a question of substitution of 0, from (37) and 
ds 


= from (33) in (32) to obtain Tj( y ,s). After rearranging 
X 








= 6S bz mb bot om ee | a 
| oe (+b)? an i ret j ic : =a = 


To invert (37) it is convenient to express Car n,s)as follows 


rm 
= “Nb -1S 
B= se ENR, (s+b2] (39) 
n=0 
where 
ee = oa 
ky = 2% BM 
kK, = oH by + Km bi b2 nH 
z (40) 
k, = Z OX mn b, Sala) 
3 
k, = a b, ban 
Similarly 
Boy Se see 
— -"b, - $ S+b -Y¥ . 
J) =se a -€ oe, (s+ bz) (41) 
h={ 





where 


ky = “ b2 
ky” ~ 2% bib2 Hh + %, ba 
k,” = bb» (Xm bi +20, + «ban ) (42) 
Re o by bz (1+2b24 ) 
Ke = ot by b. "1 
For the inversion it is useful to employ the formula 


a n=! 
£\s-e | = (£)? (a) = 


for integer n 2 1. For n=0, the formula is made valid by adding 
the term 5 (t) to the right side; this is the Dirac delta function. In 
order to express the result in terms of modified Bessel function of 
first kind order zero and one only, we use the following recursion 
formula 
TS) = Ly ls) = (ayes (44) 
for modified Bessel functions (Cf.W. G. Bickley [5], Number (20a)). 
It is also found convenient to introduce the function 
anes hy 1 (2 Yat) (45) 
for which the following differention formulas are valid. 
dh(¢ nt) _h( 24/2, n-1st) -(b+ /t+Y2t) iC , nist) 
dt (46) 
=h( @ +472, nt1;t) -(b+ °/t-9/2t) hi @, nt) 


In order to invert (39) consider the following function 
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F (48) = éyhe : F (1.t)= t"|§ ws) (47) 


using the formulas, 


£ le]. §(t) + {Z Ty (2 {at ) (48) 


and 
/ wi 
£ se]. (=) : reece) n> (49) 


yields, 


n- 


n-l 


F(n. 5) = k, §(t) ‘Yh (a) ; Tn (2494) (50) 
Next using the recursion formula (44), (50) may be expressed as 


mot) =k § (t) +m, 1(2 Yat +m, {= I, (2 {at) (57) 


where 





= t (52) 
m, = x, (/t) tk, + 2 Ke - k3) +kg(/a) 
ae iy Ale aa 


4 = ky + (Ya) (k,” esas) 


Thus we have 


~| = t -n = bit 
od |: ) bless] e° Sct) + me 15(2 (OF ) + (53) 
{= a 1, (2 {@#)| 
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Now let 
G(n,S+bz) =S-f(,3+b2) (54) 


Thus 


g(t) = a [a (nssets) -£ [s-Ferasery] (55) 
= Lt" fF (or) + f. s f ("1 .S+b2)- F(o)| 
= k(t) + 2 F (ot) 


To perform the indicated differentiation, (46) is used; thus 


Spat / 
G(n,t) =e 4080 + (2by- ba) § (E) +C, 1, (2 fat) + (56) 


Cy 1, (.{at)] 


where 
Co= bbat bnl (b-2b) a bl obit (bi-2b,+ =} (57) 
n 
ye - (* /byb2)(a/t) 7+ 9 (2b: bi) (a/e) + 
(2x bi by by -% by b"- a bboy - 2% b, [ny + 2%,b1 »)(t/a)*- 
(rq by b n) (t/a)? 


It is now easy to obtain the inversion of (39), which may be written 


as 


a6 





= -nbi -ns 


Q,(1>S) — oe ay ee eae S + b2) (58) 
Thus 
O(n, t) = O for ten 
-y bi (59) 
a,(1,t) =€ G(n,t-n) 


proceeding ina similar way, consider the function 


F as aes 
fiqis)= & \ fe S" 5 F(vted f (ws) (60) 
Thus 
n-! 


Pinot) =) (£) * Tn (2 (oe) 


Expanding the sum in (61) and using the recursion formula (44) yields 


Fiinst) =m I, (2Vat) + fE m T, (2 (ot) 


(61) 


(62) 
where 
t 2R4 6h. a yh 
Mo = hy + (£)( 43 - - = )+(4) Rs 
¥ a he ry ay A : ri (63) 
mus (K-23 y 24 88 ) 4 (E) (hae 2) 
Again let 
G™ ("55+ b2) =Sf"(1,S+ b2 ) (64) 
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Thus 


G"(n,t) = rd fe (15+) = an sf*(n, S+ 2) 


LF (us) + £°(SP*(noseh)-F(vd)] (65) 


iW 


+ 5 
t = 
1 o(t) + ) F (7 2 t)} 
The indicated differentiation in (65) is obtained using (46). Thus 


¥ -bat 
G (x,t) =€ ' 


% 
x bn S(t) +C, L,(2Vat) + CG (ovat) (66) 





where 


Co = (xb 2bib2 XN — %. bz 
0 = (%,bo + 2bibs %- % be 2) + (t/a) %mbibe | 
L 2 
° ( 26, A ee) b, - ben +2b2+ 8/n — byba % ox )- 
2 
(t/a) % bibs nN al 
a y 
Ch opto (8) [ob tbba eb] 
- 
+ X bi big (b,-2b2+1/n) (t/a) 7 


+ (% e i n} (by =) \oyt -3/n) (t/a) ” 


It is again easy to obtain the inversion of (41), which may be written as 


== -nb, —“S 
Ti(qes)= @ +e -G™ (1, S+b2) (68) 
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Thus 


for 16 AL 


{ 
° 


Ti(1> +t) 


ay 
g 4 G(q,t-n) for tz 


Ny 


Ti (n> t) 


Although the inversion of ey (nH ,s) and Dal nm ,s) were performed in 
as compact a manner as possible, it is desirable to rewrite 6,( 4 ,t) 


and Ty ( v\ at) as 


A(no.t)=Ti(yt) =o fae ae) 


-bi, -b,T 
A@(nt) = e ae ee, 8c) 4 Xn (2b,-b2) §(c) + 


Co1.(2 (at) + C11 (2127)| for tpn (aly) 


and 


-nby -b,T 


Ti(v-th=e -e og ybeSte) +601, (e/me)+cneeT) (72) 


for t2N 


It is of interest for future discussion to obtain the limiting values of 


O,(4,t) and T,(n,t)as T-+ o* 
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rior ("3 tee ‘Ln co. ‘5 NE) « (07) 


T+ ot Tr o* 


~b," “ 
ae de ; Peete bin] (abs) rent OT+...) 
—>0 


; ( ut I Ge (ar)* a 


byb2T 4 {Tv ‘ 


Thus 


, mae 
- O(n,t) = e ra baba | 3 bait rt (2bi-ba) 


4 di Km 
-2 + 73 
ie | (53) 


And similarly 


Lim Ty (not) = a nn [or (USE) + €F L215) 


To 

b. (74) 
-€ lwb ie ban 
= oe 2+ 2byba Xm -o4 b, "+ ve: bi 


Since the solutions obtained in (71) and (72) are complicated, they are 

organized so that the character of the solution can be seenas clearly 

as possible. Equation (70) shows that there is no disturbance ahead 

of the wave front. It is also noted that at the wave front (T =t-= 0) 
0, and T, are singular, both containing delta functions and @ yin , t) 

being dominated by the derivative of the delta function. Thus 9,(,t) 


2 
is not of order unity and € @, is not smallincomparisonto @ (.,t) 
2) 


at the wave front. 
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Although it is evident that the labor of obtaining higher order 
corrections (that is, higher order terms in the perturbation expansion) 
would be most formidable, it should be equally evident that such 
additional corrections would have a nature similar to that of the 
corrections already obtained; namely, they would exhibit discontinuous 
behavior at the wave front and would fail completely ahead of the wave 
front (that is, they would be consistent in predicting identically zero 
temperature). However, it is to be expected that they would possess 
calculable limiting values, from the left, as the wave front is 
approached. In effect, we thus have what appears to bea perfectly 
valid way of determining the behavior behind the wave front, but one 
which fails at, and ahead of the wave front. The nature of the failure 
of our original perturbation expansion is characteristic of singular 
perturbation problems [9], and we are therefore led to seek a singular 
perturbation expansion which will resolve the difficulties which have 


been encountered. This is done in the following subsection. 


C. SINGULAR PERTURBATION ANALYSIS 
We have seen that the perturbation solution represented by the 
series (12) and (13) involve no difficulties, except those of analytical 
complexity, forany Tyo , that is, at any point behind the wave 
front. However, this solution is useless for points ahead of the front 
™ <o , since it predicts zero temperature for fluid and solid at all 


points ahead of the wave front; physical considerations indicate that 
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diffusion through the wave front will result in non-zero temperatures 
ahead of the front. A mathematical indication of the difficulty may be 
seen from the fact that ae 6, is not small at the wave front, and 
presumably the same thing could be shown for oe O° etc... Accords 
ingly, we now look for another solution which will be valid ahead of the 


wave front. Todo this it is convenient to introduce the new variable 


en (tt) =x- Vt (75) 
and we now seek new solutions @ ( ¥ oo Pree | E; t) which will be valid 
for ¥>0 

Making the change of independent variable from (",t) to (¥ aE) 


provides the differential equations 


2 
og 7 et (8-7) oe 
eh 
Be. nog te eb (o-7) (77) 


Note that the transformation (75) changes our viewpoint so that, in 
effect, we are now moving to the right with the front so that the fluid 
appears to be stationary and the pipe appears to be moving to the left 
with velocity V. 

Note that in contradistinction with the situation dealt with pre- 


viously, now € (to the first power) appears in the governing equations. 


qo 





Accordingly € , not see appears as the natural parameter with respect 
to which we hope to construct a perturbation expansion. Thus we deal 


with 


O(E,t) =O (Et) + € (Et) A(Et)+-- (78) 


Micrtlee duieot) + 6 1 (Sst) + ¢ T, (est)+--- (79) 


(Later when we ''match" these new solutions to the solutions discussed 
in the previous subsection, we will add an asterisk to distinguish the 
new solutions. 


The conditions that all these new functions must satisfy are 


Shy ——p 0 5) pe O as CO (80) 
6-=To =0 for t=o. |, (81) 
moet = 0, 1,2, 


Additionally, they must satisfy appropriate conditions at the wave 


front. Actual selection, or determinations, of these conditions for 





$ = 0 will be made in the next subsection where the matching process 


| will be treated. 


Substituting (78) and (79) into (76) and (77) comparing powers of € , 
_ gives rise to the following sequence of equations. For the fluid, we 


will have 
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v e 
x @, = & + bi ( 0,- To) (82) 


¢ ° 

* 0, 6 + by (0, -T) (83) 
4 ° 

“6, = + b; (0, -T) (84) 


and for the solid we obtain 


[, =o (85) 
7 ¢ 
“1, =-T +1, -b2 (@-To) (86) 


«1, =-h + Ty - b2 (8-7) (87) 


(where the prime denotes partial derivative with respect to §& and 
dot denotes partial derivative with respect to time). 


Condition (80) and equation (85) show that 


Ts = 0 (88) 
Now consider the first non-trivial equation in the sequence, 


namely (82). This is a classical diffusion equation, the solution of 


which, subject to (80), is 


Bit 
O,(¥.t) =A, e » erfe (3 /2fat ) (89) 
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The constant A, is available to satisfy a boundary condition for € =0. 
The next equation to be considered is (86) which in view of (88) may be 


written as 


¢ 
Ti a ol b O; (90) 
Together with conditions (80) and (81), we may write down the solution 


of (90) as 


00 


=b)b 
T)(§,t) = Ao 2 € feck (§/2{«,t ) de (91) 
g 


Note that no new constant of integration is required in the last step. 

The next equation is (83) which is similar to (82) except that Th» 
given by (91), is involved. The solution could be obtained but we will 
not do more than remark that another constant of integration will be 
introduced. Then (87) provides I>; no new constant is introduced. 
Next (84) provides a means of determining § 3; another constant of 
integration is involved. 

Briefly the solution for 0; (iol 2 ee inviolwessasconstant 
A;, which is actually a function of time t. Since the equations for T. 
feeewtirst order, their solutions do not require the introduction of 
arbitrary constants or functions, equation (80) and (81) providing the 
necessary conditions. 


No difficulties are encountered (except those of increasing 


analytical complexity) in continuing this process and we expect 
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convergent series (78) and (79). However, as yet we have made no 
move to relate these expansions to any conditions at the left which 
might serve for the establishment of the as yet undetermined functions 


A (t). We proceed to do this in what follows immediately. 
i 


D. MATCHING OF SOLUTIONS 

We now employ an asterisk to distinguish the solutions (@ ‘ and ea 
of the immediately preceding subsection, valid ahead of the front 
(i.e., for E >o ), from the solutions ( oO; and T.) of the earlier 


subsection, valid behind the front (i.e., for €>0). We will employ 


the following notations. 


Dib pe (x,t) a Q(t) ry (92) 


Xx—>Vt- 


i 


Zan W(t) = ke {t) = Tr (93) 
x—» Vt 
-* —* 
Eco) ar acer) ae 
Fo 
— x =e 


le Va Sak (95) 


Tees) 


(Note that here an overbar does not indicate a Laplace transform. ) 
Examples are 


Oe eee (96) 
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1 Aa ome 8, (97) 


- Z 
0, =e bit ajbbrt | bb, t + (2b,- b2)t + 2-2 + “| (98) 


f 
T. = e bit (X&b,42b,b,. «Kt x bt AAS 
l 7S ( %, ae pe £ ie $ rae £ (ee (99) 
T, = 0 (100) 
x -b,t 
a Siem (101) 


@ 


EG 
= A,b, @ | enge (S/2fat ) dg (102) 


4| 
r 3 
| 


The functions 9.” and ie contain as yet undetermined functions of t, 
namely Ao = A(t), A, = A,(t), etc., and we must be able to evaluate 
these systematically. Our procedure is to evaluate so that the per- 
turbation sequences to the left and right of the front, truncated soas 
to include only terms containing € to anarbitrary even power (and 
lower powers) will match at the interface. 


The first matching is to require 


- —~ * 
@, = % (103) 


= 


To 


aa | 


(104) 


Equation (104) is automatically satisfied. Equation (103) provides a 


condition which can be used to determine the constant Ao: We find 
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fee 2X (t) = I. Thus, our first (matched) approximation is given by 
i), (31), (89), and (91) with AG = 1 in the latter two. 


The second matching is 


_ ae -* = — # 

Sree = @ 2eb o & E (105) 
= oe — — * pp . 
Thté TT = to teh + eé tz (106) 


ed 


* x sk bad 
The functions ) , Oo, : eT , and T. have introduced two additional 


unknown functions Aj(t) and A, (t). Equations (105) and (106) provide 
conditions for their determination. 
Similarly the next matching is 


2- 4 _ _ = 22 ee. 
B46 0+6€0=2=6 tea +60, +65 468, (107) 
— * —* 2-—* 3 +. 
bre h+eT, = % + €T +€ Th +€7T3 4+ € Ty (108) 


in which two additional functions, A2(t) and A, (t) appear and which 
afford means for their determination. Proceeding in this manner, all 
functions A.(t) may be determined. At each stage we have matched 
solutions to the same degree in € , on either side of the front. 

By this procedure, eventually (1) and (4) are satisfied. Moreover, 
at any stage, the necessary coefficients A,(t) are calculable so that 
the approximation affordstruncated series , available both to the left 


and to the right of the front, with matching at the front. Thus it appears 
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that we have succeeded in describing a process which provides two 
sets of series expansions, one valid behind the front and the other 
valid ahead of the front. (A '" et'' consists of a solution @ for the 
fluid and a solution T for the solid.) The labor of constructing the 
successive functions that occur in these series is considerable and 
we have actually obtained only a few of them. We have actually made 
the "first matching" and could, with additional labor, determine 
sufficient additional functions to accomplish second, and even third 
matchings. However, the writer confesses toa feeling of uncertainty 
in this entire matter of constructing the ordinary and the singular 
perturbation expansions and the matter of performing a matching 
where it is necessary to join them. Accordingly he suggests that it 
may be fruitful to consider this matter further, employing the insights 
which might be provided by a more sophisticated mathematical 


treatment. 


Pee DIFFERENT SOLUTION AHEAD OF TEMPERATURE FRONT 


The solution given in this subsection may be of some interest 


f 


although it involves a certain basic inconsistency. In this solution 

we are concerned with the temperature variation ahead of the front 
and are able to consider both radial and axial variation in the fluid, 
as well as both axial and radial diffusion in the fluid. However, we 
assume that the temperature of the solid is uniformly zero and that 


there is perfect heat transfer through the interface; that ig, we 
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consider an isothermal boundary condition for the fluid. (It would 
not alter the problem profoundly if a finite surface coefficient of heat 
transfer were assumed.) Thus, this problem is not a conjugated 
problem. For a left end boundary condition we assume that the fluid 
temperature is a constant, independent of radius, which depends on 
the location of the front. Accordingly, the degree of approximation 
and inconsistency in this solution is considerably greater than is to 
be found in other problems treated herein. However, the results may 
be useful in considering the temperature distribution ahead of the 
front at quite small values of time, before significant radial variation 
has built up in the front. 

We will consider equation (II-1) and (II-2) with isotropic thermal con- 


ductivities and subject them to a change of variable 


S=x-Vt : t=t (109) 


Then equation (II-1) reduces to 


ya 2 
08 108,98 = | 38 (110) 
ye ore m- dt 


Similarly equation (II-2) reduces to 


ip Si ee, UIE I ie 
det ae Rew" eh oe = x St (111) 


We note that the effect of the change of variable, (109), is to change 


our view point as follows; the fluid remains stationary and the pipe 


IB! 7 





wall moves to the left with the same vclocity V of the fluid. If the 
velocity V is large then we state the following assumption 
Assume that for €>0 the pipe wall is essentially at zero 
temperature. 
Hence we will concentrate our attention on equation (109) and 
find the solution when equation (109) is subjected to the following 


conditions 


O=4 ¥ =0 » oO frga , t>o fells) 
B= 0 Eo , o<rsa , to (113) 
B=0 OSE Lam ;, r= A , t20 (114) 
28 =o oS E Kw , rao » tzZo (ley 
n= 0 oS E<am , o<rla , t=o (116) 


(Note that ¥=o0 implies x=Vt, and §€-» oo impliesx—>™ ) 


Assume the solution to be 


9 = ) As > RI (Clay, (117) 


Nal 


condition (114) yields 


J,(Ana) =0 (118) 


and condition (115) is already satisfied. 
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Next, taking Laplace transforms of (110) and (117), we find 


o>, <a g aL = 
re (An + 2) ~, =0 (119) 


A solution of (119) which satisfies (113) is 


?, = B, exp (- +E) (120) 


Thus, writing C= AaB, » we have ° 


9 = Y Co Ora, exe | RE (121) 


Condition (112) then yields 


0o 
== Y Cn, (rat) (122) 


n={ 
Using the orthogonality of Bessel functions, viz. 
a 2 
(oe 
[QF (rot) Ply SE ll (An) (12 3:) 
a] . 


Meone= tr and zero if n # m, we find, 
Gas oe ae a (124) 


Hence 





= AG) ; an 
= fe exp {- 
; (0,9) J (An 2) E Cs ne) _ 





Inverting (125) we find 


me - 

Shee 2% dat Dad 

Q(E,r r,t) = Dee Sint) Caer (SS ee vente (Exitehet 26) 
(9 J, (20) se eNKse 

(Note ; that solution (126) is for the case when the wave front is 


maintained at unity.) It is now possible to use Duhamel's super- 


position integral when the condition at § = 0 is 


| _bit 
Q, (0,1, t) She 0, (it) a (ia 
To 
which is the simplest estimate we have, Cf. (96). Thus, we obtain 
t 
= bi ~ 
O(E,r,t)= eo 32 c be ie 2)| 4s (128) 
0 Eats 


where the term in brackets may be obtained by differentiating (126). 
Obviously, this solution is for cylindrical geometry, but the case 


of slab geometry differs only in detail. 
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Peer coo rt WHEREsrLUID TEMPERARUR Ease ERAGE pea) 
AXIAL CONDUCTION IN CONDUIT IS NEGLECTED 


SLAB GEOMETRY 
1, Development of Solution 
In this section equations (II-19) and (II-20), together with the 
appropriate boundary and initial conditions will be further developed. 


For convenience of reference these equations are rewritten below. 


~e Vv 26 _ | 26 ft - 

ee x gt ax, (O°) - 

¢ ; f 

oe «6 OT (2) 
dye x ot 

O=o of X <00 > t=o (3) 
O=1 X=0 y t20 (4) 
Go X00 » to (5) 
=. o<g <b » t=0 (6) 
ole u=b 5 tzZo (7) 
oy 

Ol a 
-K3Te 4 (9-7) u=0 5 tz (8) 


Note, that 9@= @ (x,t) appearing here is the average fluid temperature 


as defined in (II-18); we have omitted the asterisk uniformly throughout 





this section. Also, T = T (x, y,t) is the temperature in the duct wall; 
we have pe nitted the subscript 2 on y; in this section, the symbol y 
denotes y, in figure 1b. We also note that this analysis neglects 
the axial conduction in the duct wall. Finally, note that the surface 
heat transfer coefficient appearing above is the "equivalent" film 
resistance coefficient p discussed following (II-17); we omit writing 
the tilde in (1) and (8) above and in what follows. 

In order to solve (1) and (2), Laplace transformation is used. 


Equations (1), (3), (4) and (5) become 


xO MOF |S pes, 

Set Ox 0 “Gk, [ (T) | (9) 
—- | 

O= = X=0O (10) 
6 —>o \ =2los (11) 


Similarly equations (2), (6), (7) and (8) become 


{ 


Oo 


Ta (12) 
“Ku 5 = 4 (8-7) 0gx<0 5 y= (13) 
oT =o ofxCa » y=b (14) 





The solution of (12) is 


td 


ee -% | 

T= A@ +Be (15) 
where 

- =a (16) 


Condition (14) permits evaluating B interms of A, and we arrive 


aa be $d 09 $(2b-y) 
eC +e (17) 


T= A 


Since the terms in parenthesis are clearly independent of x, the 
imposition of condition (13) implies that A in (17) must be a function 


ofxands. Thus using (13) to express Qin terms of T, we find 
a= 7 (s)-A(x,s) (18) 


where 


ms)= fe 4 (1-28) (148) 9 


Substituting from (17) and (18) into (9) yields a homogeneous 
second order ordinary differential equation for A(x,s). After re- 


arranging and substitution for f3 (s) this becomes 


ala: ae = 4(s). A(x;5) (20 
r 





where 


s)a So LL (Km t Bank (64) 
ds) See (GUI Sexerew 7) _ 


The general solution of (20) is 
A(x,s) = cong (.) +g) i} + D.exP Vex 5) a0) (22) 
as f 


where C and D are constants of integration which depend onS. The 
boundary conditions (10) and (11) in conjunction with (18) define the 


following conditions for A (x ,s) 
A (x55) = as = 0 
Ca === 0 X —> a (24) 


Condition (24) clearly implies that D = 0 whereas (25) yields 


| 
S./3(s) 


Therefore 
Vx 
A(x,s)= —— Te exp| aa, RV A ri) (26) 


Substitution of A(x,s) from (26) above into (18) and (17) yields the 
fluid and pipe wall temperature distribution in Laplace transform form; 


V1Z. 





A Ge 1 Km % tanh (64) 
as ME yh) +, (a) 7 


T(x,435) = cosh [g(b-y)} . 
cosh (bg) + Ree (bg) Abwh (b 4) (28) 


je ls 
z Km ¥ tank ( 63) 
+} exp| TS -x (Sete (4,5) +8 (x05 
l “ ME AK Nt Ss (bg) tanh (63) = Fi ) 
In order to invert (27) and (28), we first consider @ (x, s). Denoting 


the radicand by m we have. 


2. VL Kw : 
(LJ § tatty (Sle sae) 


yar was “ (bf) Lawak ( b¢ 





mM = 


It will be convenient to cast this equation in another form. Write 


he -'% 2-i fS (30) 


so that 


S = -x ) +7 yi }\ (amp 
Then (29) becomes 


nm = Ba [(X, J S27] -ba tonnbl (x) Be -%2 4] 
B,, — br tan (bd) 





(32) 
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where 


| 
| 


Km —? (33) 


are Biot numbers for the wall and for the fluid, respectively. We 
are concerned with finding the singularities of (27). Except for the 
pole at s=0, these will all be branch points corresponding to the zeros 
and poles of m . Considering the numerator and denominator of (32) 
separately, we first note that the denominator has poles for br = 
(n+ 1/2) 41 ; n=0,1,2... . These will also be poles for the 


numerator unless 


2 


2 wan 2a 
Ee: - xe, tt Mal Ee (34) 


We may assume that this equality does not hold for any integer n 
for otherwise we could make an infinitesimally small change in any 
of the physical parameters so as to assure that the equality is violated, 
without making any significant change inthe physical problem. Thus, 
the only zeros and poles of m correspond to the zeros of the num- 
erator and denominator, respectively. 

The condition for zeros of the numerator (i.e., zeros of m ), is 


that 


tan (Ab) =f, (a) (35) 


where 
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)" an 
= (a) rs " (36) 


f(>) = b (3 KIB Be —¥ 


and the condition for zeros of the denominator (i.e., poles of m ), 


is that 
tan(bd) = f, ( d ) (Oe 
where 
f£(rX) = Sm (38) 
bA 


The real zeros of (35) and (37) may be found easily. Their location 
is illustrated in Figure 3. The relative position of the dashed line 
and the dot-dash line in Figure 3, is important. For small X itis 
evident that £ ne) < f, ( » ). Also fi I » ) has a vertical 
asymptote for a value of }\ greater than that of its single zero. Thus 
the positions are as shown in Figure 3. 

Finally, we must assume that there are no other zeros than the 
real zeros indicated in Figure 3. However, both (36) and (38) satisfy 
the hypotheses of Appendix B which shows that there are no non-real 
zeros. 

Thus it is clear that the only singularities of m are poles; 
furthermore, m_ has only real zeros and poles, infinite in number, 


and alternating with each other, of which that having the smallest 








Figure 3: 


ae 


Poles and Zeros 


(Stab Geometry) 
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= tan(bdA) 
f, 0A) 
#04) 


x = poles 


of m(s) 





numerical value is a zero. If we designate the abcissae of the inter- 
sections denoted by 0 in Figure 3 by the symbol Ne and those 
denoted by x by the symbol we ,» we conclude that the singularities 
of (27) are: 
a. branch points (corresponding to zeros of m ) at 
SVG ee \e and 
m “An : 
b. branch points (corresponding to poles of m ) at 
Sey2 
S=- x » is 
c. a pole at S=0 with residue l 
Furthermore, the branch points of type a and type b alternate with 
each other, that closest to the origin being of type a. 
Thus the hypotheses of the inversion theorem given in Appendix A 
hereof are satisfied. Accordingly, we may immediately write the 


inverse as follows 


OF int 
Sle) nile see) [é Ne kp (u)] an (39) 
where B(x) = (\m(-n)| ) 


Similarly T(x,y,s) as given by (28) is conveniently organized in the 
form of the product of two Laplace transforms. Accordingly, its 


inverse may be easily obtained in the form of a convolution integral. 
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Tope 4) Sede ah) cs fa bey 


2d,,b + Ain(2 yb) 
=] t 


-X Nt fx >? (40) 
2 fe Ta(xendn 


°o 


where An are the roots of (37) and the integrand contains the 
function (39). We have used the formula 


-| } S 
cosh (as%) 


CosR (bs) + & (bS14) Anh (65%) 
(41) 


00 Baer 
mM _ Pu Din (dn) Cos (Anb) 


2,6 + Aw (2bdn) 
ri=| 


which may be derived directly or inferred from equations given in 
IZ} examples 5-3 and 7-11. 
2. Asymptotic Approximations 
Thus the desired solution of the system described by equation 
(1) through (8) has been achieved. We will not in this work discuss 


computational aspects, but will now obtain some asymptotic 


approximations. 


Equation (27) can also be written 


= 1 ve. Vv Km 
8(x,s) = pla 2A, X ee “aan X/p, +X /ph 2) (42) 


120 





42 2,4 3,6 a8 
ae > 55sec 


Consider the radicand to be expanded in ascending powers of s. 


2 
O(x,s) = + exp m =x (x eae (ee OG) (44) 


%¢ 
Considering small values of s, which corresponds to large t, 
the first asymptotic approximation is obtained by neglecting all powers 


of s except the zeroth power. Correspondingly 


6ix,5) = (45) 


Le 
S 
which yields 


§(x,t) = 1 (46) 


This is evidently correct as an asymptotic expansion since for any 
finite x the temperature approaches unity for sufficiently large time. 
A nontrivial asymptotic expansion is obtained by retaining 


the first power of s. Correspondingly 


a) - J Vx : o Kern $ 


The reason for writing A, /A¢ , the area ratio, in place of 
b/a which naturally would appear in the last term of the radicand will 


appear in Section VB2. 


12] 


The inversion of (47) is tabulated [29], thus 


vn 
Q(x,t) = = erfc (S48) + Of erte Ge ed (48) 


where 


V"= v/( +fe5. A /p-5 4, ) (49) 


is a retarded velocity which is identical to that obtained by Munk [24], 


% 
> + [RA] BEA, (50 


is an effective thermal diffusivity for the fluid. This solution should 


and 


be compared with that given in (VII-10) where a different set of 
idealizations leads to a similar result. There, the retarded velocity 
(VII-5) is the sarne as (49), but the diffusivity (VII-6) is slightly 
different from that given in (50). 

For the corresponding solution for the conduit wall temp- 
erature we employ (48) in (40). 

Retaining the second power of s would yield a useful 
improved asymptotic approximation if one could find a simple explicit 
inversion formula. Unfortunately nothing simpler than the general 
inversion given in Appendix A has come to light, and this has provided 


the complete solution given in (39) and (40). 


Ea 


Be CYLINDRICAL GEOMETRY 
1. Development of Solution 
A similar treatment to that given in Section VA is employed. 
Only the significantly different details will be treated explicitly. 


Equation (2) will be replaced by 


WT, twee! wT 


{ 
Ge r or Km Ot el 


with the corresponding conditions 


T =o ; asér¢b ) t=o (52) 

ol - r=b »  tp2o (53) 

or 

-K, ot ak (o-T) r=a ; t> oO (54) 
r 


replacing (6), (7), and (8). Remarks similar to those following (8) 
are applicable here. 
As previously, Laplace transformations are made with 


respect to t. The solution of the transformed version of (51) is 
T = AL, (Gr) + BKo(4r) (55) 
where 


f =  S/% (56) 


is the same as previously. Here I. and K, denote modified Bessel 


functions. 


E23 





Proceeding as in VAI, one arrives at equations similar to 
(27) and (28) except that the function tan (qb) appearing therein is 


to be replaced by another function which we denote 


T, (x) Ki (84) = 1, (8x) Ki (x) 
T, (8x) Ki(x) + Ko(Bx) Ty (x) 


(ip, 


bank (x, 6) = 


(The initial b denotes ''Bessel.'' The remainder (anh) is intended to 
indicate a similarity to the function tanf(x) .) 


Corresponding to (29), we have in this case, the equation. 


2 2. 


= 2% : be 8G) 
m Q 
i ce (63) bank (b2 , /b) 
Corresponding to the well known relation 
ys tank (x) Ee ae (2x) (59) 
it is not difficult to show that 
dL boak(x,@) = bam (4x, 8) (60) 
where 
Te (x)Y, (8%) Ye (x) (8x) 


ban (x,4) = 
( T(x) Yo (0%) -Y, (x) J, (4x) 
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Thus corresponding to (32), upon defining 


Na - 44 


we obtain 


ie 
velo cas a 
Bm - (b») bon (bd, a/b) 


where 


_ bh a 
Ba J a= 


(Note the factor 2 in the present definition of Be ) 


(62) 


(63) 


(64) 


From this point, the development is precisely the same as in VAI. 


Considering the counterparts of equations (35) and (36), precisely the 


same treatment is applicable. Figure 4, shows a graph similar to 


Figure 3. Figure 4 was actually constructed by computer analysis in 


order to show that the function ban (x, 3 ) behaves like tan(x); the 


value of @ used in constructing Figure 4 was fi = U.9.. 


There is one difference which we would prefer to ignore. 


In the previous development, the fact that m had no non-real zeros 


or poles was proved by the argument contained in Appendix B. To 


construct a similar proof for the case in which tan x is replaced by 


ban (x, #? ) appears terribly formidable and we simply omit doing so. 


Thus, a solution corresponding to (39) and expressed in 


exactly the same notation as (39) (that is, identical to (39)) is 


R25 








\ 
~ | \ 
i 
| | 5 
eos SS 
| oe } + 
| | i + |! | | 
J 
N a | ee 
= Si = 2 by 
| 
\ | 
N 
\ 
\ 
P| 
\ | 
\ 
—_ 
\ 
—— y = ban(bi,f) 
----y = f, 0) 
—-e es Y = #50) 
oO - zeros, X - poles 


Figure 4: Poles and Zeros of ms) 


(Cylindrical Geometry) 
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applicable in the present case. However, the function F (r,s) in (28) 
must now be replaced by 

E (r,s) = Pe aR 

[- Seg Liga) + Le (59) lk (44) + [He g Ki ga)+K, (59)]1, Gi) 


The inversion of (65) may be obtained using a solution discussed 


5 ) 


in Section 13.4 of [7], by first specializing to the case of a unit step 
input, and then by differentiating the result with respect to time. The 


result, corresponding to (40), is 


— matt 
Tirx,t) =k we LE, (6aspn) + “E06; | case” fx ‘ean a 
dA, wom er ion au ae (x; 
Ti(%,0) * £/}? 


where pn are the roots of 


ep = a6 (a,b; p) 
ho G (ab ¥) 





and 


E49) = F(X) Yn(04) - Ye (6%) Tn (24) (68) 


&¢. Asymptotic Approximations 
It is possible to proceed here in exactly the same way as in 
Section VA2. The principal difference is that in the present case of 
cylindrical geometry, in place of (43) we now define xX differently 


and arrive at different coefficients in the expansion. Thus 


ae 





X = (bq) bank(be a/b) = = (P- a)S +... (69) 


where... denotes terms in 5%, etc. It requires painstaking work 
to establish that the coefficient of s° vanishes and that the first 
power term is as in (69) but the calculations are not illuminating and 
are not given here. 

As in the case for slab geometry, the first non-trivial 
asymptotic expansion is obtained by retaining the first power of s 
in (69); correspondingly, we arrive at equation (41) the inverse of 
which is given in (48). In writing (47) we anticipated the present 
development by writing A, /A, in the final term of the radicand soas 
to be able to accommodate to both slab geometry and cylindrical 


geometry with one equation. 
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VI. .A, GENERALIZATION OF A PROBLEM GONSIDE Ra. 
BY Ree 


This case is described by (II-10) and (II-14) subject to conditions 


that will be stated below. The problem description is thus 


xe x, = 3 em — ™ 
a erat (9-T) (2) 
Q@=] X =0 » t2o (3) 
@=0 X—»ya~ , t>30 (4) 
@=0 o<x <a 3 t=o (5) 
T=0 0x <a 5 t=o (6) 
6 = 9(x,t), T = T(x,t) are average temperatures in the sense of 


(II-7), (II-8), but we omit using the asterisk which distinguished this 
significance. The heat transfer coefficient h in (1) and (2) above is 
defined by (II-9). 

This problem is a generalization of that considered in Section 
7-5 of [2], the generalization coming about by the inclusion here of 
axial conductivity in the fluid, a mechanism not considered in Arpaci's 


analysis. This problem is a specialization of the problem considered 
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in Section IV hereof since here we neglect axial conduction in the 


duct wall, whereas this mechanism is included in the analysis given 


in Section IV. 


Taking Laplace transformation with respect to time t, one obtains 


i hag pay _— 
OND N08 AS aah, 1 = 

. x, sr (z+ ‘) 8 + bi] =0 (7) 
ST = b2(6-T) (8) 


where 





Combining (7) and (8), one obtains 


26 V0 _/s.,4_ bbz ye _ 10 
aye x. Ox zt arr, ) 8 = 0 (10) 


3 a = X=0 (11) 
=s6 Nee aes (12) 
Thus, the solution of (10) may be written 
is eee oe ier ere an 
Ee L exp] YE x yr Sa by- bbe LS) 
> 2x. 2%, ay S+ b2 


an increasing exponential having been rejected in order to satisfy 


(12). We write (13) as 
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bre 8 Vx xX  /(S4r) (St) 
at 2K, a (S+ bz ) eee 


: f 





where 
\= ets v 22.) (15) 
f 
2 
‘odie Un a oo 8 sas 
1 
(17) 


V 


W= b2+b, & re 
2 wee 2 a 


It is easy to show that A> p Se However, we want to indicate 


and la relative to by . To see this, consider 


the magnitudes of xX 


the function 
2 2 2 
S V b2 bz V 
f(s) = oX +*S( gat bh +e) + 4X2 se) 


f 
the vanishing of which corresponds to the roots p Ne Ue viz 
We note that the coefficient of the quadratic 


Gee) = f(-p ) = 0. 


and we note that 


term is positive, 
f{(-b>) = -byb2 <0. (19) 
Thus, we easily see that, 
(20) 


A>br>p>o 


Thus the condition of the theorem stated in Appendix A are satisfied 


and we can write the inverse of (14) 


13] 





bz 
—Kt 


O(x.t)=1- 7 exp (YS ) nr p (ny}an (21) 


, “st 
-[£ him (x p(n} 


PN 
Fe (1) = H(\-p )(A-11) a (22) 
%, (2-4) 
Using the convolution theorem to invert (8), we obtain 


t 
rks 
T(ut)=bz/e To (x, tn) dy (23) 


© 


where 


We should remark that equations (1) and (2) above may be 
obtained by specializing equations (61) and (62) of Hiep's treatment 
[13]. However, we note again that Hiep's results depend upon 
numerical inversion of Laplace transformations, a process which 


Oureexperience indicates is of doubtful accuracy and reliability. 


2 


VII. AN IMPROVEMENT ON A RESULT BY W. MUNK 


Under the assumption that radial conduction is so efficient that 
essentially the same temperature prevails at all points of a particular 
fluid and metal cross section, W. Munk [24] obtained a very simple 
approximate solution for what he called ''The Delayed Hot Water 
Problem.'' In the present section we obtain a very simple improve- 
ment over Munk's result. Hiep [13] has also obtained the improvement 
given here and he mentions that it has been obtained by others. Our 
reason for presenting it is that it is a natural deduction to be drawn 
from developments already made in this thesis. 


We repeat equation (II-13) and the corresponding conditions. 





ey (1+ Kou Amy BT y ot = (I+ Km Am %e ) QT 


Pe Ae x2 OK VT KA aq Ot e 
T=1 at x=0 3» t20 (2) 
ee) as X—»00 9 t 20 (3) 
Tee at OfX<O » +t2=0 (4) 


We have used T to denote the common value of average temperature. 


By defining 


ae (fy eau 
V ev [(ie Ber) (5) 


x ae P, Bono 
XH = (1 a) X¢ / (+ Be) (6) 
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we may write 


2, x 
oe Ly, Oly eel 


7 
ox? axe OG ae 
Taking Laplace transforms with respect to time, we obtain 
*aA XAT = | 
of oT _y OW Soar (8) 
Ox OX 
The solution of (8) satisfying all conditions is 
* 
= ¥ \2 
20 V x V S 
T == exp axe axe) + KH 
Inverting (9) by the use of standard tables [29] we obtain 
*% 
vt LES % 
2 2x t 
T (x,t) = ih erfc (a=V Ee) + € erge( X40 | (10) 
2 2 fart 2 yxkt 
This is the generalization described in the title of this section. 
Munk's result is 
T (x,t) = H (x-yt) (11) 
maere HH ( ) denotes the Heaviside unit step function. Thus 


the present result (10) is like Munk's in predicting a temperature 
wave front moving with a reduced velocity We but, unlike Munk's 
' it also provides for a softening of the wave front due to axial 


diffusion in fluid and pipe. 





VIII. DISCUSSION 


Pe COPE OF ANALYSIS 

The general problem as formulated by equations (II-1) and (II-2) 
for the case of cylindrical geometry and equations (II-3) and (II-4) 
for the case of slab geometry, even with the simplifying assumptions 
of subsection [IA2, are mathematically very involved. These 
equations contain two dependent variables, namely the fluid and the 
pipe wall temperature each of which is a function of three independent 
variables. There is also a formidable number of boundary and initial 
conditions to be satisfied. In Section III the problem is reduced to 
the solution of an integral equation the only unknown in which is the 
temperature at the interface, which is a function of only two variables; 
no auxiliary conditions are required. However, as may clearly be 
observed by the reader, the integral equation of interest is math- 
ematically very involved, and retains a great share of the complexity 
of the original problem, for the following reasons. 

Firstly, the unknown function appears only under the integral 
sign (the so called integral equation of the first type). Secondly, 
there are two integral signs under which the function of interest 
appears. Thirdly, one of the integrals has a semi-infinite range, 
namely o to oo. Fourthly, the kernels of integrations are 


complicated functions. 
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The complexity of the general problem, either in its original 
form or as compressed into this integral equation is such as to have 
led us to focus attention upon further simplifications which are more 
tractable from a mathematical viewpoint and the analyses of which 
are treated in other sections of this thesis. However, the integral 
equation formulation is believed to be of significance and its treat- 
ment, probably by numerical procedure, is considered feasible. 

Further simplifying assumptions are introduced in IID, in order 
to reduce the mathematical difficulties. Ina search for the various 
methods which would render the general problem amenable to math- 
ematical analysis, it is found that by defining a suitable mean temp- 
erature in the direction perpendicular to the direction of flow, and 
by making a further assumption of constant surface heat transfer 
coefficient, the useful formulation given in (IJ-10) and (II-11) is 
obtained. 

Assumption of large velocity justifies the use of a perturbation 
technique, which is presented in Section IV hereof. 

The zeroth order terms, which corresponds toa solution given 
in Section 15-3 of [7] and also in example 7.5 of [2],are obtained, and 
first corrections to the fluid and pipe temperatures are obtained 
explicitly. These solutions are found not to be valid at and beyond 
the temperature front. Accordingly a singular perturbation tech- 
nique is used to obtain an expression which traces the behavior ahead 


of the front. 
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In Section Vl a further physical assumption, namely neglecting 
axial conduction in the solid, leads to a problem for which we are 
able to find the closed form solution by using the method of Laplace 
transformation. This affords an improvement over the solutions given 
in [7] and [2], cited in the preceding paragraph. 

In Section V the physical assumptions that are introduced consist 
of considering an average fluid temperature and neglecting axial con- 
duction in the solid. This case is solved using the method of Laplace 
transformation. Even though the methods of transform calculus 
reduce the partial differential equations to a system which can be 
readily solved in the transform domain, sometimes delicate math- 
ematical handling is required to invert the solutions back to the real 
variable domain. Sucha case is that of equation (V-27) and (V-28), 
where the solutions in transformed form contain an infinite number 
of singularities (branh points). An elegant theorem for this inversion, 
established by Professor J. E. Brock, is presented in Appendix A. 
Thus closed form solutions are obtained. Furthermore, use is made 
of the properties of Laplace transform to yield some useful asymptotic 
solutions. 

In Section VII the view taken is that of W. Munk [24]. It was found 
possible to improve upon his solution by accounting for the axial 
conduction in the fluid as well as the solid wall. This treatment 
assumes, as Munk did, an instantaneous common constant uniform 


temperature at any cross section of the fluid and wall. Although 
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this result is obtained independently, it turns out to be known; see 
Appendix II of [13]. This solution is useful for small fluid velocities; 
Cf-Hiep [13], p. 32 and Preston [27]. 

We will make no further attempt to relate the problems considered 
herein to practical engineering situations. Such references as [7] 
and [2] are devoted to problems having similar relations to practical 
engineering situations. Their importance lies in the fact that they do 
provide definite mathematical solutions to well posed mathematical 
problems which have relation to the physical world of heat transfer 
engineering; it is the obligation of the person who employs such 
results to investigate the degree to which the mathematical conditions 


are satisfied in his own physical application. 


B. RECOMMENDATIONS 

Most important is treatment of the integral equation developed 
in Section III. Remarks about this are included in Sepeccmon JG bleh, 

It would be useful to develop asymptotic formulas to employ in 
place of the exact inversion in cases where there are severe numerical 
difficulties associated with the latter. The method of steepest descents 
(e.g. Cf. Carrier Crook and Pearson - Section 6-3 [8]), appears 
promising for this purpose. 

There are several mathematical details which should be filled 
in (although in treatments similar to our own they are usually glossed 


over). These all relate to locations of the zeros and poles of 
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transcendental functions. In the case of equations (V-59) and (II- 155), 
we have simply assumed that the zeros and poles are all real; under 
this assumption they are not particularly difficult to locate. How- 
ever a proof that there are no complex zeros or poles is really 
meourred , 

Inasmuch as instantaneously higher rates of heat transfer are to 
be expected at and near a rapidly moving front, it might be e296 Clea 
that an analysis employing finite speed of heat propagation would be 
enlightening. Accordingly, it is suggested that such an analysis be 
attempted in one of the simpler subcases dealt with herein; Cf., for 
examples, Section 1-6 of [10]. 

Also since numerical difficulties are to be expected in applying 
the inversion theorem of Appendix A to the heat transfer problems 
herein, it is suggested that actual numerical calculations be under- 
taken. A multiple-precision facility, such as described in [14] and 
[28] seems necessary for this purpose. 

Another thing which comes to mind is to relax the specialization 
herein because of the assumption of slug flow, i.e., constant fluid 
velocity. It seems clear that the assumption of velocity as a given 
function of the radial or transverse dimension would provide definite 
differential systems for analysis, but that these would, in general, 
be much more difficult to deal with than those considered herein. 

The treatment given by Hiep [13] relies heavily on the accuracy 


and validity of a numerical inversion of Laplace transform due to 
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Salzer [30], [31], [32]. However our experimentation with such methods 
of inversion showed that reliably accurate results could be obtained 
only in limited ranges and for relatively simple functions. This 
persuaded us to look for other analytical methods of treating the 
problem, such as perturbation method treated in Section IV hereof. 
However, it is probable that further development of techniques of 
numerical inversion may provide a useful and powerful tool for 
analyses such as those studied herein. 

We should also mention our feeling that no great difficulty would 
be associated with inferring response to a unit step change in left 
end temperature but retaining the other conditions specified in the 
analysis of Sparrow and DeFarias [33], who employed a sinusoidally 
varying left end boundary condition. 

Also although throughout we have assumed that the conduit had 
perfect exterior insulation so that there was no heat flux through this 
surface, it seems clear that other exterior boundary conditions (such 
as an isothermal condition) would, in most applications considered 
herein, involve only changes in detail and not essential changes in 


procedure. 
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APPENDIX A 


INVERSION OF SOME LAPLACE TRANS FORMS 


The theorem stated and proved in this appendix provides the 
Laplace inversion of (V-27), (V-28), and (VI-14) of the text. This 
appendix is largely a paraphrase of [6]. 


THEOREM: If the function to be inverted is of the form 


f(s) = = exe (ymca) | ae) 


where m(s) is a meromorphic function whose zeros s = -a; and whose 


poles s = -b, satisfy 
ON <b Qn bp --- Ca fb, +++ o@ (2) 
=| 
Then, for t > 0, the inverse Laplace transform Fes f fF¢s) 


is given by 


N ar 
oetmop a fapoln 
K=| an 


where 


B(n) =/P" Gn) (4) 


The interpretation of the sum in (3) is as follows. If there is a 
finite number N of pairs a, ¢ b,, no special interpretation is 


necessary. If there is a finite number (N-1I) of pairs a g DT and an 
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felaitional zero s = - ant, an >> bN_-}> an Nth pole s =~- by = 
is to be added soas to make N pairs. If there is an infinite number 
of pairs a, < bes take N = © 

In proving the theorem stated above we consider a sequence {R | 
of positive numbers, R, # ay b, for any nand any k, R_—-+o as 
n—»o00, and we consider a corresponding sequence of closed contours 
co in the complex s- plane each of which encloses the pole of f(s) at 


s = 0 but excludes all of the branch points at s = -a, ands = - by. 


k 
The contour C., consists of four subcontours Cy}, Cn2, Cn3 and C,y, 
as follows (see Figure A-1l). Cn} is a vertical line froms = - iR_: 
tos = iR,, indented as shown by a small semicircle passing to the 
16 

memo: s =~ 0. Cp isthearcs =R,€@ , Teo ct and Ci4 is the 

, 2 = = 

es) 16 . 
similar arc s = R,€ DAUSTBS . Finally C,3 is the contour which 
passes to the right froms = -Rn, along the negative real s-axis to 
S$ =- a} with small semicircular indentations into the upper half s - 
plane about each zero (-a,) and pole (-b;) and then passes to the left 


from s =- a, tos=- ism with small semicircular indentations into the 


lower half s- plane, as indicated in Figure A-1l. 





Figure A-1l: Contour o- 
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Clearly by Cauchy's theorem 


st ) 
T= ap Ge fis) ds = eel Vee) | one re) 5) 


~ 2re 
Cr 


forn=1,2,3... . Also, by the complex inversion theorem for 


Laplace transforms 
ii 1 | é Nels [5 (6) 
— a emmeatl —— a 
Cat 


as Nn — »00 * 
To proceed further we wish to establish an inequality for arg 
[m(s)]. To do this we consider complex vectors emanating from the 


points s=-a, and s="by, and we write 


{+2 =Uep(ig.) > Ueto ; <> eT (7) 
Qk = ak 
1+ . =u, exP(ty,) | Udo 3 SH <T (8) 
K 
im) m) 
Yn (s) =m(o) [| (hte jae (i+ =) (9) 


= p(s) exp (i Pn) | (10) 


where 
(n) 


(v1) 
pe (s) = ple) Tl U+ TT UR (11) 


and (n) 


Bn = ) % Vy (12) 


The sums and products written with a superscript (n) indicate that 
only those terms are to be considered for whicha, < R, or on < Ro» 


as may be appropriate. 
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Now note for any s in the (closed) upper half plane we have 
% 2 ae 3 Sas yea! ees (13) 


A geometrical argument is given in Figure A-2. Now either 


there are N pairs a,» b, each less than R» so that 


a=) (4%) 2 a4 
K=/ 
or there is an extra, unpaired zero, so that 
N-I 
n=) ( 4-%) + P, >o (15) 
K=| 


In the first case we can write 


Bn 


(16) 


1 
ae, 

! 
s 
= 
<< 


=~b, =- Qa, 


Figure A-2: Demonstration that f. > Vv, 
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whereas in the second case we can write 


By = fo vy, 2h) o>? <Tt (17) 
Ks/ 
Thus in either case 
oS fp. < TT (18) 


Similarly by considering an s in the lower half plane we find 
WS ffs Se (19) 


so that in all cases 


-1S An. ST (20) 
Thus, if we write 
st 
IT = I ds 
Siler 7a ie f(s) (21) 
Cn2 
we note that 
eS) || eels tare 22 
| | oe Pn (22) 











exp |-Pa os(Pr)] gb 


\- 


uniformly on C,92, by virtue of the inequality (20). Thus, by Jordan's 
Lemma I,2 —»0O as n-—~»oo. Similarly the integral on C,4 vanishes 
as n—>»oo 

Next we consider the integrals around semicircles surrounding 


poles and zeros on contours C,3 The result is obvious near a zero 


st -a,t 
s=-a; @ we » | Role since m(s) a , and 
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the length of the arc goes to zero as its radius does. Thus the integral 
on such semicircles about a zero, s = - ay vanishes with the radius. 


However, for the integral on a semicircle about a pole, s = - by, we 


must invoke (20) again. We have 


exe[-Fms] = &xp | Pu Cos ( Par)) <! (23) 





even though 4 i—? 0 as the radius of the indentation goes to zero. 
Thus the integral on such a semicircle about a pole, s = - br, vanishes 
with the radius. 

We remark that although we have reason to believe, on physical 
grounds, the correctness of an evaluation by Kantola [18], we feel 
that his argument intended to show vanishing of such an integral ona 
semicircle surrounding a branch point is deficient. Our argument, 
based on (20), is not applicable to his case. 

If now we write (5) in the form 

Tui = 4,7 'n2 - In3 - 14 (24) 
and let n——» oo _ , and also let the radii of all semicircular indenta- 
tions on C,,3 go to zero, we find 

F(t) = exp [- p (0) [Po (25) 


where 





ig & | = 2x p (st - Ym(sy) ds (26) 


21re 
C3 
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and C, now is a contour extending from s = - 00 along the upper edge 


3 
ofa branch cut, to the right, passing around s = - ay and then 
extending along the lower edge of the branch cut tos =- © . The 


branch cut is the negative real s - axis from - © to-a,. Thereare 
no longer any semicircular identations in Cz. 

To reduce (26) to a more convenient form, we first introduce the 
function V (x), for real x, which denotes the excess of the number 
of zeros (s = - a,.) over the number of poles (s = - by) to the right of 
the ant s =x. It is clear that Y=0 in the interval -a,4] << s 4 -b, 
while vV=1 in the interval “by gs < =a, Writing (3 now, rather 
than n , since we have taken n—>o and all zeros and poles 
are included in our accounting, we note that on the upper edge of the 
cut fp = VT! /2 while on the lower edge 3 == yif2. We also new 


write p rather than Pn . Thus 


-Q, 


— fs enp [st —p exp (-vma)] ds + 


"07 
= 00 


ft exp| st _p exp (-i vm) ds 


“2 


(27) 


-a 
i] st 


= [ee levee ou] ct mr giom 


if » = 0, the quantity in braces vanishes, while if v= 1 it 


becomes 


exp(-tH) -exp(dp) =- 24 Atm p (28) 
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Thus 


wea Palme 


Kat bx 
If we now change the variable of integration from s to "N=- s, we 
obtain 
N “nh 
ar haet aeee B (nl dn (30 
T/ | 
kK=\ Ak 
and, thus, finally 
bx 
fet 
F(t) = exe | po] - i) | nin PUa) dn (31) 
K=| Ar 


The use of the upper limit of summation is consistent with the remarks 


following the statement of the theorem, the proof of which is now 


complete. 


In equation (VI-14) the solution obtained in Laplace transform 


form is a special case in which we have 
2 2 
ens) = k (s+2i) (sea) / (5 +b) (32) 


where e@ > 0:and o <Q, <b, <Q . The numbers aj, az and 


b, are known explicitly by (VI-9, 15, 16). We have 


H(o) = Po = &\/& Ar /bi (33) 
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and the inversion is 


by 


be ae 
me) = € - —- - 
Q 
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APPENDIXGE 


NONEXISTENCE OF COMPLEX ROOTS OF EQUATION (V-36) 


Demonstration that such equations as (V-36) and (V-38) have only 
real roots are occasionally called for in order to employ Cauchy's 
Theorem (and Jordan's Lemma) to evaluate the inverse Laplace 
transforms encountered in heat conduction problems. Carslaw and 
Jaeger [7] have some general remarks on this matter (Section 12.8 
IV , page 324), making use of Sturmian theory; their reference is to 
Chapter X of [15]. Polya and Szego [26] consider similar problems, 
most of them concerned with the zeros of polynomials, but some,e.g., 
Auf.172, of nature similar to ours. They use various methods of 
proof. Professor J. E. Brock has constructed the following elementary 
proof for the case encountered in equation (V-36). The case encountered 
in equation (V-61) wherein the simple trignometric tangent function is 
replaced by a generalization involving quotients of expressions involv- 
ing J- and Y- Bessel functions, is obviously much more difficult and 
we simply make no attempt to construct a proof inthis case, taking 
refuge in the (weak) argument that the function we encounter is 
"similar'' to the tangent function. 


Thus, we consider 


C(n-§ - x") 
n- x* 


xtanx = C>o 5, n>o of8Kn (1) 


Let x = a+ib, a, b real and presume b # 0 


We will show that this leads to a contradiction 


tan (a + ib) = (tana +itanhb)/ ( 1 - i tana tanhb). 


Denoting the left side of (1) by L, then 
b 


lta tan x — (a + 7b} —_____—___—,, | | —__________ 
\- tana tank b]\ 145 Tan a lank 


which may be written as 


_ (AAnra-b Adhib)+é (bAm2a +a acd 2b ) ee 


x tan x 
D 


where 


D. = cos2a + cosh2db. 


] 
Next considering the right hand side of equation (1) and denoting it 


by R we obtain 
c(t 6-a8)-2606] (nas B) pacad 
(nab) 2606 (nate 

(mt b* at)’, 4a°b-s (nb a4 )-24 ab § 
: 2. 


R 


where 


D2 =[(nebh a?) +4a°b |/< 


At this point it is well to observe that for b #0, coshb > 1 


While for any reala, cosa 2 - 1 


and thus we conclude that 


Di > 0 


(n-a*s b? ) + ae 
(4) 


(5) 


(6) 





Also 





2 * 
Modo y= Ltip 
Des = ( s > 0 n>o  , b>. (Za) 
nowe(n - a~+b2) and (Zab) cannot both vanish. 

Now if we equate the real and imaginary parts of (3) and (4), 


we obtain 


. 2 
Win ra-bAmk2b (ue Ween Kh lnnle a) 





= (8) 
D, D2 
by equating the real parts, and 
Dont Noi 2 b 
Ae (9) 
Dy 2 
by equating the imaginary parts. From (9) above we may write 
. z 
Apn2a __ Q? Panh2 b 26a (10) 
Dy b Dy D2 
and substituting in (8) we obtain 
2 2,2 2 2 Z 
DAS alka) 4a b-§(n+b-2 )416a 


b Di D2 


Now considering the numerator on the right of (11). Callit N, and 


consider it as a function of § - then 


a =-(n+b-2) + 20° <30 <heeeee (12) 


| gay 


First, supposing that fe > O, then the minimum value of N occurs 


for = 0. 
a 2,2 
Nee (nab - a ) +4ab > 0 
Next suppose that iF < 0, then the minimum value of N occurs for 


oe n. 


7s 
Noni = (n+ esas) =f 4a bx (n+ aa) tan q* 


a ,2 é 24 ee) 
=n(@eb)+(aeb) > o 
Thus inany event MW > 0, and since for any real b # 0 we have 
Joon ty 2 
Domb 2b > Ors (14) 


b 
thus the left side of (11) is less than zero whereas the right side is 
greater than zero; which is a contradiction. 
Thus we cannot have roots of the form x = atib witha and b real, 


and b #0. Therefore we can only have real roots. 


v5 





uO), 


a. 
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Sidered under the assumption of constant uniform velocity and constant material 
properties. The general problem, including all significant heat transfer mech- 
anisms, except radiation, is reduced to an integral equation containing only one 
unknown function, namely an interface temperature as a function of time and | 
Mxial position. By neglecting one or more heat transfer mechanisms, several 
Other formulations are obtained. These are analyzed so as to provide approximate 
and closed form solutions not previously available. Similar formulations and 
results are obtained for the case of fluid between parallel plates. 
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